시계열 분석을 위한 1번째 단계이다. 이번 단계에서는 독립 변수가 하나인 Univariate 시계열 분석( ETS( Error, Trend, Seasonal ) 모델과 같은 AR, MA, ARIMA )에 대한 여러 주제를 살펴보자.
모든 회사는 심한 경쟁, 연구 실패, 노사 관계, 인플레이션, 침체, 정부 정책 변화와 같은 여러 Challenge와 위험을 마주한다. 당연하게도, 모든 기업은 불확실성과 위험에 노출되어 있다. 위험으로 인한 피해를 줄이고, 다가오는 위험은 어떤 것인가를 알기 위해 예측이 필요하다. 다양한 예측 모델이 있는데, 대부분 많이 사용하는 것은 아래와 같다.
왜 예측마다 다른 테크닉을 사용할까? 데이터의 생김새, 특징이 제 각각 다르기 때문에 서로 다른 테크닉을 사용한다. 예를들어 회귀나 CART는 여러 개의 예측 변수( predictors )에서 1개의 응답( response )를 얻는다.
예측은 historical data를 이용해서 미래 Trend의 방향을 결정하는 예측 가능한 정보를 만들어 내기 위한 기술이다. 경영에서는 미래의 예산과 지출을 어떻게 산정할 지 예측할 때 사용한다.
이번 단계에서, 시계열 데이터의 상세 내용과 어떻게 분석, 예측할 수 있는지 알아봅시다.
시계열이란 시간의 흐름에 따라서 모은 동일한 변수들의 측정값( 그냥 measurement라 하겠다. ) 급수( 역시 그냥 Series 라 하겠다... )이다. 이러한 measurement들은 일정한 시간 간격에 따라서 만들어진다. 시계열은 시간의 순서에 따른 Data point들의 Series이다. 대부분, 시계열은 해당 구간 내에서 동일한 point 간격으로 떨어진 형태로 보여진다. 따라서 discrete-time data의 sequence 이다.
시계열 데이터의 간격
시계열 분석은 자산, 증권, 경제 변수가 시간에 따라서 어떻게 변했는지를 보여주는 유용한 방법이다. 선택한 데이터를 동일한 구간에서 다른 변수와 비교하고 얼마나 변화가 생겼는지 확인할 때도 사용된다.
대부분 우리가 마주치는 시계열 데이터의 대표적인 사례는 아래와 같다.
다시 말해서, 우리는 시계열 데이터를 수집하고, 데이터 기반으로 의사 결정을 내리는데 활용할 수 있다.
다음에 언급된 특징들은 시계열 분석 및 머신러닝 분석을 쉽지 않게 만든다.
데이터가 cross-sectional인 경우, observation의 순서는 중요하지 않지만, 시계열의 경우에는 중요하다.
시계열 데이터는 몇 가지 공통괸 추정치를 갖는다. 가장 유념해야 할 것은...
아주 장기 예측은 잘 동작하지 않는다. !!
사용가능한 데이터가 있다고 해도, 우리는 몇 주기 이상의 예측을 시도해서는 안된다.
import numpy as np, pandas as pd, matplotlib.pyplot as plt, os
DATA_ROOT = os.path.join( os.getcwd(), 'TimeSeriesTutorial' )
for dirname, _, filenames in os.walk( DATA_ROOT ):
for filename in filenames:
print(os.path.join(dirname, filename))
C:\Users\user\Jupyter\TimeSeriesTutorial\air_passengers\AirPassengers.csv C:\Users\user\Jupyter\TimeSeriesTutorial\shampoo_sales_dataset\shampoo_sales.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\AirPassengers.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\AirPax.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\AirTemp.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\BOE-XUDLERD.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\Champagne.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\CrudeOil.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\daily-min-temperatures.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\daily-total-female-births.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\Emission.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\EngWage.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\Gasoline.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\GDPIndia.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\GDPUS.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\MaunaLoa-1.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\MaunaLoa.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\Petrol.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\pollution.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\PortugalPort.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\RetailFood.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\RetailTurnover.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\shampoo.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\TractorSales.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\wage_growth.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\WaterConsumption.csv C:\Users\user\Jupyter\TimeSeriesTutorial\time-series-data\WeeklyClosing.csv
### 단일 변수 시계열 데이터 읽고, 상위 몇 줄만 확인해 보자.
univariate_series = pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'daily-min-temperatures.csv' ), header = 0, index_col = 0, parse_dates = True, squeeze = True)
univariate_series.head()
Date 1981-01-01 20.7 1981-01-02 17.9 1981-01-03 18.8 1981-01-04 14.6 1981-01-05 15.8 Name: Temp, dtype: float64
# 패턴좀 확인해 보기 위해서 plot 좀 해 보자.
univariate_series.plot()
plt.ylabel('Minimum Temp')
plt.title('Minimum temperature in Southern Hemisphere \n from 1981 to 1990')
plt.show()
보통 우리는 다음 몇 달치 온도를 예측하려고 할 때, 과거의 값들을 보고, 어떤 패턴이 있는지 찾아내려고 한다. 매년 어떤 계절적 특성 패턴을 살펴보자. 특정 관찰값들은 미래 값들을 예측하는데 도움을 준다. 주의점: 우리는 "온도"라는 단일 변수( 과거 19년치 )만 사용한다는 점을 유의하라.
이런 분석을 Univariate Time Series Analysis/Forecasting( 단일변수 시계열 분석, 예측 ) 이라고 한다.
다중변수 시계열 데이터는 한 시점에 1개 초과의 데이터를 가진다. 각 변수는 과거값 뿐만 아니라, 다른 변수의 영향도 받는다. 이런 종속관계가 미래의 값을 예측하는데 사용된다. 아래 웹 사이트를 한변 살펴보아라.
대기 오염 예측
이 데이터셋은 지난 5년간 중국 북경 주재 미국 대사관에서 매 시간마다 측정한 날씨, 대기 오염도 데이터이다. 본 데이터는 시간, PM2.5, Dew point( 이슬점 ), 온도, 기압, 풍향, 풍속, 누적 강우량, 누적 적설량을 포함하고 있다.
Raw 데이터 형태( 전체 Feature )는 아래와 같다.
| Sl No | Variable | Description |
|---|---|---|
| 1 | No | row number |
| 2 | qyear | year of data in this row |
| 3 | month | month of data in this row |
| 4 | day | day of data in this row |
| 5 | hour | hour of data in this row |
| 6 | pm2.5 | PM2.5 concentration |
| 7 | DEWP | Dew Point |
| 8 | TEMP | Temperature |
| 9 | PRES | Pressure |
| 10 | cbwd | Combined wind direction |
| 11 | Iws | Cumulated wind speed |
| 12 | Is | Cumulated hours of snow |
| 13 | Ir | Cumulated hours of rain |
이전 시간의 날씨, 오염 상태가 있을 때, 다음 시간의 대기 오염 수치에 대해서 예측해보자.
첫 24시간 중 PM2.5 데이터 중 NA 값이 있는지 확인해보자. 첫 번째 Row를 제거해야 할 것이다. 그 외, 몇 가지 NA 값이 있을 것이다. 우선 그 NA 데이터들을 0으로 채워보자.
# 날짜 정보 처리
from datetime import datetime
parse = lambda x: datetime.strptime( x, '%Y %m %d %H' )
# 데이터 로드
pollution_df = pd.read_csv(
os.path.join( DATA_ROOT, 'time-series-data', 'pollution.csv' ),
parse_dates = [ [ 'year', 'month', 'day', 'hour' ] ],
index_col = 0, date_parser = parse
)
pollution_df.drop( 'No', axis = 1, inplace = True )
# 수동으로 컬럼 이름 변경
pollution_df.columns = [ 'pollution', 'dew', 'temp', 'press', 'wnd_dir', 'wnd_spd', 'snow', 'rain' ]
pollution_df.index.name = 'date'
# NA 전부 0으로 변경
pollution_df[ 'pollution' ].fillna( 0, inplace = True )
# 첫 24시간 데이터 삭제
pollution_df = pollution_df[ 24: ]
# 첫 5개 데이터 확인
print( pollution_df.head( 5 ) )
pollution dew temp press wnd_dir wnd_spd snow rain date 2010-01-02 00:00:00 129.0 -16 -4.0 1020.0 SE 1.79 0 0 2010-01-02 01:00:00 148.0 -15 -4.0 1020.0 SE 2.68 0 0 2010-01-02 02:00:00 159.0 -11 -5.0 1021.0 SE 3.57 0 0 2010-01-02 03:00:00 181.0 -7 -5.0 1022.0 SE 5.36 1 0 2010-01-02 04:00:00 138.0 -7 -5.0 1022.0 SE 6.25 2 0
values = pollution_df.values
# 컬럼을 plot, Wind speed는 ordinal 데이터가 아니라서 제외
groups = [ 0, 1, 2, 3, 5, 6, 7 ]
i = 1
plt.figure()
for group in groups:
plt.subplot( len( groups ), 1, i )
plt.plot( values[ :, group ] )
plt.title( pollution_df.columns[ group ], y = 0.5, loc = 'right' )
i += 1
plt.show()
airPax_df = pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'AirPassengers.csv' ) )
# str -> datetime으로 파싱
airPax_df[ 'Month' ] = pd.to_datetime(
airPax_df[ 'Month' ], infer_datetime_format = True
)
airPax_df_indexed = airPax_df.set_index( [ 'Month' ] ) # "월"을 Index로 설정
airPax_df_indexed.head(5) # 확인
| #Passengers | |
|---|---|
| Month | |
| 1949-01-01 | 112 |
| 1949-02-01 | 118 |
| 1949-03-01 | 132 |
| 1949-04-01 | 129 |
| 1949-05-01 | 121 |
plt.plot(airPax_df_indexed)
plt.show()
### Dataframe을 ts1.csv 파일로 저장
airPax_df_indexed.to_csv('ts1.csv', index = True, sep = ',')
### 다시 읽어서 제대로 불러들였는지 확인
series1 = pd.read_csv( 'ts1.csv', header = 0 )
print( type( series1 ) )
print( series1.head( 2 ).T )
<class 'pandas.core.frame.DataFrame'>
0 1
Month 1949-01-01 1949-02-01
#Passengers 112 118
인도 GDP Series 데이터를 읽고, TS 데이터를 저장해 본다. 데이터는 1960-01-01 ~ 2017-12-31 까지이다. 아래는 Python으로 객체 데이터 Binary로 IO하는 Pickle에 대한 설명 자료 이다. 참고.
Python의 어떤 객체라고 Pickle을 통해서 IO가 가능하다. Pickle은 Python 객체( list, dict 등... )를 character stream으로 변환하는 방법 중 하나이다. 본 방법을 이용하면 다른 Python script에서 동일한 형태의 Python 객체를 불러올때 유용하게 사용할 수 있다.
TS 객체를 저장하고, Pickle object를 불러와보자.
india_gdp_df = pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'GDPIndia.csv' ) )
date_rng = pd.date_range( start = '1/1/1960', end = '31/12/2017', freq = 'A' )
india_gdp_df[ 'TimeIndex' ] = pd.DataFrame( date_rng, columns = [ 'Year' ] )
india_gdp_df.head( 5 ).T
| 0 | 1 | 2 | 3 | 4 | |
|---|---|---|---|---|---|
| Year | 1960 | 1961 | 1962 | 1963 | 1964 |
| GDPpercapita | 81.284764 | 84.426437 | 88.914919 | 100.048592 | 114.315161 |
| TimeIndex | 1960-12-31 00:00:00 | 1961-12-31 00:00:00 | 1962-12-31 00:00:00 | 1963-12-31 00:00:00 | 1964-12-31 00:00:00 |
plt.plot( india_gdp_df.TimeIndex, india_gdp_df.GDPpercapita )
plt.legend( loc = 'best' )
plt.show()
No handles with labels found to put in legend.
### Pickle로 데이터 IO 테스트
import pickle
with open( 'GDPIndia.obj', 'wb' ) as fp: # India GDP dataframe 데이터를 binary 형태의 pickle로 저장
pickle.dump( india_gdp_df, fp )
with open( 'GDPIndia.obj', 'rb' ) as fp:
india_gdp1_df = pickle.load( fp )
india_gdp1_df.head( 5 ).T
| 0 | 1 | 2 | 3 | 4 | |
|---|---|---|---|---|---|
| Year | 1960 | 1961 | 1962 | 1963 | 1964 |
| GDPpercapita | 81.284764 | 84.426437 | 88.914919 | 100.048592 | 114.315161 |
| TimeIndex | 1960-12-31 00:00:00 | 1961-12-31 00:00:00 | 1962-12-31 00:00:00 | 1963-12-31 00:00:00 | 1964-12-31 00:00:00 |
시계열의 구성요소는 아래와 같다. :-
추세( Trend ) :- 장기간에 걸친 상대적 상위, 하위 값들의 점진적인 변동 및 움직임
1. 점점 올라가는 것으로 보이면, 상승 경향( Uptrend )이라고 한다.
2. 점점 내려가는 것으로 보이면, 하강 경향( Downtrend )이라고 한다.
3. 추세가 없으면, 수평( horizontal ) 혹은 변동 없음( Stationary ) 경향이라고 한다.
계절성( Seasonality ) :- 상승, 하강의 Swing 상태를 의미. 정해진 시간동안 어떤 추세가 반복된다. 예를 들어, 춥고, 더운 날씨에 산다면, 에어컨 이용 가격은 여름에 높고, 겨울에 낮게 나올 것이다.
샴푸 판매 데이터셋을 예제로 들어보자.
본 데이터셋은 3년 동안의 월별 샴푸 판매량 데이터를 나타낸 것이다. 판매량이 단위이고, 36개 관찰결과가 있다. 원본 데이터셋은 Makridakis, Wheelwright, and Hyndman (1998) 이 만들었다.
Data source:https://github.com/jbrownlee/Datasets
아래가 첫 5개의 데이터 이다.
| Month | Sales |
|---|---|
| 1-01 | 266.0 |
| 1-02 | 145.9 |
| 1-03 | 183.1 |
| 1-04 | 119.3 |
| 1-05 | 180.3 |
parser = lambda x: datetime.strptime( '190' + x, '%Y-%m' )
series = pd.read_csv(
os.path.join( DATA_ROOT, 'time-series-data', 'shampoo.csv' ),
header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser
)
series.plot()
plt.show()
일간 최저 기온 데이터셋
본 데이터셋은 호주, 멜버른의 1981 - 1990에 걸친 10년 간의 최저 온도 값을 갖고 있다.
단위는 섭씨이며 3,650개의 측정 결과가 있다. 본 데이터의 출처는 호주 기상청이다. 5개의 데이터만 보면 아래와 같다.
| Date | Temperature |
|---|---|
| 1981-01-01 | 20.7 |
| 1981-01-02 | 17.9 |
| 1981-01-03 | 18.8 |
| 1981-01-04 | 14.6 |
| 1981-01-05 | 15.8 |
series = pd.read_csv(
os.path.join( DATA_ROOT, 'time-series-data', 'daily-min-temperatures.csv' ),
header = 0, index_col = 0, parse_dates = True, squeeze = True
)
series.plot()
plt.ylabel('Minimum Temp')
plt.title('Minimum temperature in Southern Hemisphere \n From 1981 to 1990')
plt.show()
이제 우리는 매월 변화량을 Box Plot으로 그릴 수 있게 되었다. 계절적인 요소가 여름부터 겨울까지 Swing으로 보여진다.
months = pd.DataFrame()
one_year = series['1990']
groups = one_year.groupby(pd.Grouper(freq='M'))
months = pd.concat([pd.DataFrame(x[1].values) for x in groups], axis=1)
months = pd.DataFrame(months)
months.columns = range(1,13)
months.boxplot()
plt.show()
이 그림은 최저 온도의 변화가 보여지는데, 남반구의 여름이 1월부터 시작되고, 겨울이 매년 중간, 그리고 막판에 다시 여름이 되어 가는 모습을 보여준다.
최저 온도 데이터셋을 연도별로 그룹핑해본다. Box와 Whisker plot이 연도별로 배치되서 상호 비교하기가 쉬워진다.
groups = series.groupby(pd.Grouper(freq='A'))
years = pd.DataFrame()
for name, group in groups:
years[name.year] = group.values
years.boxplot()
plt.show()
연도별 큰 차이는 보이지 않는다.
트랙터 판매 데이터 본 데이터는 2020년 매월 판매된 트랙터의 개수에 대한 데이터이다.
tractor_df = pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'TractorSales.csv' ) )
tractor_df.head(5)
| Month-Year | Number of Tractor Sold | |
|---|---|---|
| 0 | 3-Jan | 141 |
| 1 | 3-Feb | 157 |
| 2 | 3-Mar | 185 |
| 3 | 3-Apr | 199 |
| 4 | 3-May | 203 |
import calendar
dates = pd.date_range(start='2003-01-01', freq='MS', periods=len(tractor_df))
tractor_df['Month'] = dates.month
tractor_df['Month'] = tractor_df['Month'].apply(lambda x: calendar.month_abbr[x])
tractor_df['Year'] = dates.year
#Tractor.drop(['Month-Year'], axis=1, inplace=True)
tractor_df.rename(columns={'Number of Tractor Sold':'Tractor-Sales'}, inplace=True)
tractor_df = tractor_df[['Month', 'Year', 'Tractor-Sales']]
tractor_df.set_index(dates, inplace=True)
tractor_df = tractor_df[['Tractor-Sales']]
tractor_df.head(5)
| Tractor-Sales | |
|---|---|
| 2003-01-01 | 141 |
| 2003-02-01 | 157 |
| 2003-03-01 | 185 |
| 2003-04-01 | 199 |
| 2003-05-01 | 203 |
tractor_df.plot()
plt.ylabel('Tractor Sales')
plt.title("Tractor Sales from 2003 to 2014")
plt.show()
위 그래프를 보면 확연한 계절적인 요소가 보인다. 2011년 기준으로 월간 변화량을 boxplot으로 그려보자.
months = pd.DataFrame()
one_year = tractor_df.loc['2011']
groups = one_year.groupby(pd.Grouper(freq='M'))
months = pd.concat([pd.DataFrame(x[1].values) for x in groups], axis=1)
months = pd.DataFrame(months)
months.columns = range(1,13)
months.boxplot()
plt.show()
매년 5 ~ 8월 사이에 높은 판매량이 일어나는 계절적인 요인을 볼 수 있다.
tractor_df.loc['2003']
| Tractor-Sales | |
|---|---|
| 2003-01-01 | 141 |
| 2003-02-01 | 157 |
| 2003-03-01 | 185 |
| 2003-04-01 | 199 |
| 2003-05-01 | 203 |
| 2003-06-01 | 189 |
| 2003-07-01 | 207 |
| 2003-08-01 | 207 |
| 2003-09-01 | 171 |
| 2003-10-01 | 150 |
| 2003-11-01 | 138 |
| 2003-12-01 | 165 |
샴페인, 음식데이터로 동일한 분석을 수행하면 계절적 요인이 있다는 것을 볼 수 있다.
통계 보델 라이브러리들은 이에 대한 Naive, Classical한 분해 함수를 제공한다. ( seasonal_decompose() ) 여기서 Additive 인지 Multiplicative인지 설정해야 한다.
seasonal_decompose() 함수는 Decomposition 해서 4개 조각을 담은 배열 객체를 반환한다.
소매 회전율 데이터에 대한 Additive model decomposition
turnover_df= pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'RetailTurnover.csv' ) )
date_rng = pd.date_range(start='1/7/1982', end='31/3/1992', freq='Q')
turnover_df['TimeIndex'] = pd.DataFrame(date_rng, columns=['Quarter'])
turnover_df.head()
plt.plot(turnover_df.TimeIndex, turnover_df.Turnover)
plt.legend(loc='best')
plt.show()
No handles with labels found to put in legend.
from statsmodels.tsa.seasonal import seasonal_decompose
import statsmodels.api as sm
decompTurnover_df = sm.tsa.seasonal_decompose(turnover_df.Turnover, model="additive", freq=4)
decompTurnover_df.plot()
plt.show()
<ipython-input-44-bec827531ff2>:3: FutureWarning: the 'freq'' keyword is deprecated, use 'period' instead decompTurnover_df = sm.tsa.seasonal_decompose(turnover_df.Turnover, model="additive", freq=4)
위 코드를 돌리면 분해 작업을 수행한다. 4개의 Series 결과가 바로 그것이다. 어떤 추세와 계절성이 확실히 분리되어 보이는 것이 보인다.
trend = decompTurnover_df.trend
seasonal = decompTurnover_df.seasonal
residual = decompTurnover_df.resid
print(trend.head(12))
print(seasonal.head(12))
print(residual.head(12))
0 NaN 1 NaN 2 13692.5375 3 13674.1875 4 13716.7375 5 13748.1500 6 13789.3500 7 13827.8875 8 13866.6500 9 13944.9125 10 14062.1375 11 14228.4125 Name: trend, dtype: float64 0 -524.613498 1 -510.713672 2 1894.960113 3 -859.632943 4 -524.613498 5 -510.713672 6 1894.960113 7 -859.632943 8 -524.613498 9 -510.713672 10 1894.960113 11 -859.632943 Name: seasonal, dtype: float64 0 NaN 1 NaN 2 -188.697613 3 149.645443 4 -58.624002 5 34.263672 6 -88.010113 7 49.745443 8 67.263498 9 -129.998828 10 -83.197613 11 -2.279557 Name: resid, dtype: float64
# 샴페인 데이터를 이용해서 Additive model decomposition을 사용해보자.
# date_rng = pd.date_range(start='1/1/1964', end='30/9/1972', freq='M')
# date_rng
# Champ['TimeIndex'] = pd.DataFrame(date_rng, columns=['Month'])
국제 항공 여객 데이터를 이용해서 Multiplicative model decomposition을 수행해보자.
airPax_df = pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'AirPax.csv' ) )
print(airPax_df.head())
date_rng = pd.date_range(start='1/1/1949', end='31/12/1960', freq='M')
print(date_rng)
airPax_df['TimeIndex'] = pd.DataFrame(date_rng, columns=['Month'])
print(airPax_df.head())
Year Month Passenger
0 1949 Jan 112
1 1949 Feb 118
2 1949 Mar 132
3 1949 Apr 129
4 1949 May 121
DatetimeIndex(['1949-01-31', '1949-02-28', '1949-03-31', '1949-04-30',
'1949-05-31', '1949-06-30', '1949-07-31', '1949-08-31',
'1949-09-30', '1949-10-31',
...
'1960-03-31', '1960-04-30', '1960-05-31', '1960-06-30',
'1960-07-31', '1960-08-31', '1960-09-30', '1960-10-31',
'1960-11-30', '1960-12-31'],
dtype='datetime64[ns]', length=144, freq='M')
Year Month Passenger TimeIndex
0 1949 Jan 112 1949-01-31
1 1949 Feb 118 1949-02-28
2 1949 Mar 132 1949-03-31
3 1949 Apr 129 1949-04-30
4 1949 May 121 1949-05-31
decompAirPax = sm.tsa.seasonal_decompose(airPax_df.Passenger, model="multiplicative", freq=12)
decompAirPax.plot()
plt.show()
<ipython-input-46-b3e79231424f>:1: FutureWarning: the 'freq'' keyword is deprecated, use 'period' instead decompAirPax = sm.tsa.seasonal_decompose(airPax_df.Passenger, model="multiplicative", freq=12)
위 코드를 실행하면 4개으 결과가 나온다. 추세와 계절성이 확연히 분리된 것이 보인다.
seasonal = decompAirPax.seasonal
seasonal.head(4)
0 0.910230 1 0.883625 2 1.007366 3 0.975906 Name: seasonal, dtype: float64
quarterly_turnover = pd.pivot_table(turnover_df, values = "Turnover", columns = "Quarter", index = "Year")
quarterly_turnover
| Quarter | Q1 | Q2 | Q3 | Q4 |
|---|---|---|---|---|
| Year | ||||
| 1982 | NaN | NaN | 13423.2 | 13128.8 |
| 1983 | 15398.8 | 12964.2 | 13133.5 | 13271.7 |
| 1984 | 15596.3 | 13018.0 | 13409.3 | 13304.2 |
| 1985 | 15873.9 | 13366.5 | 13998.6 | 14045.1 |
| 1986 | 16650.3 | 13598.4 | 14183.2 | 14128.5 |
| 1987 | 16380.7 | 13512.8 | 14022.1 | 14231.8 |
| 1988 | 16737.0 | 14004.5 | 14165.5 | 14203.9 |
| 1989 | 16895.1 | 14248.2 | 14719.5 | 14855.8 |
| 1990 | 17361.6 | 14585.2 | 14873.5 | 14798.4 |
| 1991 | 17115.2 | 14284.9 | 14558.8 | 14914.3 |
| 1992 | 17342.3 | NaN | NaN | NaN |
quarterly_turnover.plot()
plt.show()
quarterly_turnover.boxplot()
plt.show()
petrol_df = pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'Petrol.csv' ) )
petrol_df.head()
date_rng = pd.date_range(start='1/1/2001', end='30/9/2013', freq='Q')
#date_rng
petrol_df['TimeIndex'] = pd.DataFrame(date_rng, columns=['Quarter'])
print(petrol_df.head())
plt.plot(petrol_df.TimeIndex, petrol_df.Consumption)
plt.legend(loc='best')
plt.show()
No handles with labels found to put in legend.
Year Quarter Consumption TimeIndex 0 2001 Q1 14.978 2001-03-31 1 2001 Q2 11.099 2001-06-30 2 2001 Q3 10.057 2001-09-30 3 2001 Q4 10.454 2001-12-31 4 2002 Q1 9.295 2002-03-31
Moving Average Smoothing은 간단하고, 효율적인 시계열 예측 방법이다. Smoothing은 시간의 흐름에 따른 미세한 변동성을 없애주는 기법이다.
Moving average를 계산하려면 본래 시계열 데이터에 Raw 측정값의 평균이라는 새로운 Series를 만들어야한다. 이동 평균은 Window size를 정해줘야 한다. Moving average 산출에 필요한 Raw 관찰값들의 개수이다.
Pandas의 rolling() 함수는 각 관찰값들을 하나의 Window에 넣어준다. Window 크기를 넣으면 이동하는 Window가 생성이 된다. 이동하는 Window가 생겨나면 Dataset이 변동될 때 마다의 평균값을 뽑아낼 수 있다.
airTemp_df = pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'AirTemp.csv' ) )
date_rng = pd.date_range(start='1/1/1920', end='31/12/1939', freq='M')
airTemp_df['TimeIndex'] = pd.DataFrame(date_rng, columns=['Month'])
airTemp_df.head()
| Year | Month | AvgTemp | TimeIndex | |
|---|---|---|---|---|
| 0 | 1920 | Jan | 40.6 | 1920-01-31 |
| 1 | 1920 | Feb | 40.8 | 1920-02-29 |
| 2 | 1920 | Mar | 44.4 | 1920-03-31 |
| 3 | 1920 | Apr | 46.7 | 1920-04-30 |
| 4 | 1920 | May | 54.1 | 1920-05-31 |
plt.plot(airTemp_df.TimeIndex, airTemp_df.AvgTemp)
plt.legend(loc='best')
plt.show()
No handles with labels found to put in legend.
temp_avg = airTemp_df.copy()
temp_avg['avg_forecast'] = airTemp_df['AvgTemp'].mean()
plt.figure(figsize=(12,8))
plt.plot(airTemp_df['AvgTemp'], label='Data')
plt.plot(temp_avg['avg_forecast'], label='Average Forecast')
plt.legend(loc='best')
plt.show()
mvg_avg = airTemp_df.copy()
mvg_avg['moving_avg_forecast'] = airTemp_df['AvgTemp'].rolling(12).mean()
plt.plot(airTemp_df['AvgTemp'], label='Average Temperature')
plt.plot(mvg_avg['moving_avg_forecast'], label='Moving Average Forecast')
plt.legend(loc='best')
plt.show()
USGDP_df = pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'GDPIndia.csv' ), header=0)
print(USGDP_df.head())
date_rng = pd.date_range(start='1/1/1929', end='31/12/1991', freq='A')
print(date_rng)
USGDP_df['TimeIndex'] = pd.DataFrame(date_rng, columns=['Year'])
plt.plot(USGDP_df.TimeIndex, USGDP_df.GDPpercapita)
plt.legend(loc='best')
plt.show()
No handles with labels found to put in legend.
Year GDPpercapita
0 1960 81.284764
1 1961 84.426437
2 1962 88.914919
3 1963 100.048592
4 1964 114.315161
DatetimeIndex(['1929-12-31', '1930-12-31', '1931-12-31', '1932-12-31',
'1933-12-31', '1934-12-31', '1935-12-31', '1936-12-31',
'1937-12-31', '1938-12-31', '1939-12-31', '1940-12-31',
'1941-12-31', '1942-12-31', '1943-12-31', '1944-12-31',
'1945-12-31', '1946-12-31', '1947-12-31', '1948-12-31',
'1949-12-31', '1950-12-31', '1951-12-31', '1952-12-31',
'1953-12-31', '1954-12-31', '1955-12-31', '1956-12-31',
'1957-12-31', '1958-12-31', '1959-12-31', '1960-12-31',
'1961-12-31', '1962-12-31', '1963-12-31', '1964-12-31',
'1965-12-31', '1966-12-31', '1967-12-31', '1968-12-31',
'1969-12-31', '1970-12-31', '1971-12-31', '1972-12-31',
'1973-12-31', '1974-12-31', '1975-12-31', '1976-12-31',
'1977-12-31', '1978-12-31', '1979-12-31', '1980-12-31',
'1981-12-31', '1982-12-31', '1983-12-31', '1984-12-31',
'1985-12-31', '1986-12-31', '1987-12-31', '1988-12-31',
'1989-12-31', '1990-12-31', '1991-12-31'],
dtype='datetime64[ns]', freq='A-DEC')
mvg_avg_USGDP = USGDP_df.copy()
mvg_avg_USGDP['moving_avg_forecast'] = USGDP_df['GDPpercapita'].rolling(5).mean()
plt.plot(USGDP_df['GDPpercapita'], label='US GDP')
plt.plot(mvg_avg_USGDP['moving_avg_forecast'], label='US GDP MA(5)')
plt.legend(loc='best')
plt.show()
이동 평균선이 원본 데이터에 아주 근접하게 붙어서 이동하는 것을 볼 수 있다.
IndiaGDP_df = pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'GDPIndia.csv' ), header=0)
date_rng = pd.date_range(start='1/1/1960', end='31/12/2017', freq='A')
IndiaGDP_df['TimeIndex'] = pd.DataFrame(date_rng, columns=['Year'])
print(IndiaGDP_df.head())
plt.plot(IndiaGDP_df.TimeIndex, IndiaGDP_df.GDPpercapita)
plt.legend(loc='best')
plt.show()
No handles with labels found to put in legend.
Year GDPpercapita TimeIndex 0 1960 81.284764 1960-12-31 1 1961 84.426437 1961-12-31 2 1962 88.914919 1962-12-31 3 1963 100.048592 1963-12-31 4 1964 114.315161 1964-12-31
mvg_avg_IndiaGDP = IndiaGDP_df.copy()
mvg_avg_IndiaGDP['moving_avg_forecast'] = IndiaGDP_df['GDPpercapita'].rolling(3).mean()
plt.plot(IndiaGDP_df['GDPpercapita'], label='India GDP per Capita')
plt.plot(mvg_avg_IndiaGDP['moving_avg_forecast'], label='India GDP/Capita MA(3)')
plt.legend(loc='best')
plt.show()
이동 평균선이 원본 데이터에 붙어서 이동하는 것을 알 수 있다.
간단하게 누락 데이터를 넣을 수 있다. dtype을 고려해서 넣는다.
import pandas as pd
import numpy as np
s = pd.Series([1,2,3,4,5,6])
s.loc[4] = np.NaN
print(s)
0 1.0 1 2.0 2 3.0 3 4.0 4 NaN 5 6.0 dtype: float64
Descriptive statistics와 computational statistical methods는 누락 데이터를 다룰 수 있도록 작성되어 있다.
예제:
시계열은 시간의 순서대로 나열된 데이터 포인터 Series이다. 대부분 사람들은 누락된 값이 있으면, 그 부근에 있는 데이터로 대체하려는 움직임을 보인다. 보통 이미 있는 데이터를 이용해서 Imputation( 결측값 대체 ), Interpolation( 보간 ) 방법을 이용하고는 한다.
누락 값을 대체하기 위한 Imputation을 하는 방법은 아래와 같다.
| Method | When suitable |
|---|---|
| Take average of the nearest neighbours( 가장 근접한 것들의 평균 산출 ) | Data has no seasonality( 계절성이 없을 때 ) |
| Take average of the seasons from two or all available years( 2개 이상의 연도에서 계절별 평균 산출 ) | Data has seasonality( 계절성이 있을 때 ) |
| Interpolate function of pandas( pandas의 interpolation 함수 사용 ) | |
| Linear interpolation( 선형 보간법 ) | Relationship in the interval of two samples is a first order polynomial( 2개 Sample 사이에서 1차 다항식이 나오는 경우 ) |
| Polynomial such as Quadratic or Cubic interpolation( Quadratic, Cubic 보간법과 같이 다항식이 나오는 경우 ) | Second or third order polynomial describes the interval between two samples( 2개 이상의 Sample들 사이에서 2, 3차 다항식이 나오는 경우 ) |
| Spline | Handles non-uniform spacing of samples( Sample들 사이에서 정형화된 무언가가 없을 때 ) |
이번 예제는 아래 URL을 참고해보자.
def handle_missing_values():
print()
print( format( 'How to deal with missing values in a Timeseries in Python', '*^82' ) )
# 날짜 만들고
time_index = pd.date_range( '28/03/2017', periods=5, freq='M' )
# DataFrame 만들고, Index 지정
df = pd.DataFrame(index=time_index);
print(df)
# Feature를 만든다. 누락 데이터를 넣는다.( Gap을 만든다. )
df['Sales'] = [1.0,2.0,np.nan,np.nan,5.0];
print(); print(df)
# 누락 데이터를 보간한다. 1번 방법
df1= df.interpolate();
print('\n1. Interpolate'); print(df1)
# Forward fill( 이전 데이터를 이용 ) 방법을 사용한다. 2번 방법
df2 = df.ffill();
print('\n2. ffill' ); print(df2)
# Backward fill( 다음 데이터를 이용 ) 방법을 사용한다. 3번 방법
df3 = df.bfill();
print( '\n3. bfill' ); print(df3)
# 보간법 사용, 연속되는 누락 데이터 중 1개만 처리, FFILL 사용. 4번 방법
df4 = df.interpolate(limit=1, limit_direction='forward');
print( '\n4. Interpolate, limit 1, ffill' ); print(df4)
# 보간법 사용, 연속되는 누락 데이터 중 2개만 처리, BFILL 사용. 5번 방법
df5 = df.interpolate(limit=2, limit_direction='forward');
print( '\n5. Interpolate, limit 2, bfill' ); print(df5)
handle_missing_values()
************How to deal with missing values in a Timeseries in Python*************
Empty DataFrame
Columns: []
Index: [2017-03-31 00:00:00, 2017-04-30 00:00:00, 2017-05-31 00:00:00, 2017-06-30 00:00:00, 2017-07-31 00:00:00]
Sales
2017-03-31 1.0
2017-04-30 2.0
2017-05-31 NaN
2017-06-30 NaN
2017-07-31 5.0
1. Interpolate
Sales
2017-03-31 1.0
2017-04-30 2.0
2017-05-31 3.0
2017-06-30 4.0
2017-07-31 5.0
2. ffill
Sales
2017-03-31 1.0
2017-04-30 2.0
2017-05-31 2.0
2017-06-30 2.0
2017-07-31 5.0
3. bfill
Sales
2017-03-31 1.0
2017-04-30 2.0
2017-05-31 5.0
2017-06-30 5.0
2017-07-31 5.0
4. Interpolate, limit 1, ffill
Sales
2017-03-31 1.0
2017-04-30 2.0
2017-05-31 3.0
2017-06-30 NaN
2017-07-31 5.0
5. Interpolate, limit 2, bfill
Sales
2017-03-31 1.0
2017-04-30 2.0
2017-05-31 3.0
2017-06-30 4.0
2017-07-31 5.0
예제 2 :- 물 사용 데이터
waterConsumption_df=pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'WaterConsumption.csv' ) )
waterConsumption_df.head()
| Date | reference | target | |
|---|---|---|---|
| 0 | 15-01-2010 | 12.0 | 12.0 |
| 1 | 15-02-2010 | 18.0 | 18.0 |
| 2 | 15-03-2010 | 22.0 | 22.0 |
| 3 | 15-04-2010 | 26.0 | 26.0 |
| 4 | 15-05-2010 | 31.0 | NaN |
# 일자 Column의 데이터 타입을 DateTime 형으로 변환
waterConsumption_df.Date = pd.to_datetime(waterConsumption_df.Date, format='%d-%m-%Y')
waterConsumption_df = waterConsumption_df.set_index('Date')
waterConsumption_df.head()
| reference | target | |
|---|---|---|
| Date | ||
| 2010-01-15 | 12.0 | 12.0 |
| 2010-02-15 | 18.0 | 18.0 |
| 2010-03-15 | 22.0 | 22.0 |
| 2010-04-15 | 26.0 | 26.0 |
| 2010-05-15 | 31.0 | NaN |
# 차트를 그리기 위해서, 누락 데이터를 포함한 컬럼을 하나 만든다.
waterConsumption_df = waterConsumption_df.assign(missing= np.nan)
waterConsumption_df.missing[waterConsumption_df.target.isna()] = waterConsumption_df.reference
waterConsumption_df.info()
<class 'pandas.core.frame.DataFrame'> DatetimeIndex: 96 entries, 2010-01-15 to 2017-12-15 Data columns (total 3 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 reference 96 non-null float64 1 target 75 non-null float64 2 missing 21 non-null float64 dtypes: float64(3) memory usage: 3.0 KB
waterConsumption_df.plot(style=['k--', 'bo-', 'r*'],figsize=(20, 10))
<AxesSubplot:xlabel='Date'>
# mean, median을 이용해서 채워보자.
# DataFrame에 새로운 컬름을 하나 만든다.
# df['NewCol']=0, 요 방법 말고, df = df.assign(NewCol=default_value) 이 방법 사용한다. ( Warning 방지를 위함. )
waterConsumption_df = waterConsumption_df.assign(FillMean=waterConsumption_df.target.fillna(waterConsumption_df.target.mean()))
waterConsumption_df = waterConsumption_df.assign(FillMedian=waterConsumption_df.target.fillna(waterConsumption_df.target.median()))
# Rolling average를 이용해서 Impute
waterConsumption_df = waterConsumption_df.assign(RollingMean=waterConsumption_df.target.fillna(waterConsumption_df.target.rolling(24,min_periods=1,).mean()))
# Rolling median을 이용해서 Impute
waterConsumption_df = waterConsumption_df.assign(RollingMedian=waterConsumption_df.target.fillna(waterConsumption_df.target.rolling(24,min_periods=1,).median()))
# 서로 다른 방법을 이용해서 보간을 수행해서 Impute
waterConsumption_df = waterConsumption_df.assign(InterpolateLinear=waterConsumption_df.target.interpolate(method='linear'))
waterConsumption_df = waterConsumption_df.assign(InterpolateTime=waterConsumption_df.target.interpolate(method='time'))
waterConsumption_df = waterConsumption_df.assign(InterpolateQuadratic=waterConsumption_df.target.interpolate(method='quadratic'))
waterConsumption_df = waterConsumption_df.assign(InterpolateCubic=waterConsumption_df.target.interpolate(method='cubic'))
waterConsumption_df = waterConsumption_df.assign(InterpolateSLinear=waterConsumption_df.target.interpolate(method='slinear'))
waterConsumption_df = waterConsumption_df.assign(InterpolateAkima=waterConsumption_df.target.interpolate(method='akima'))
waterConsumption_df = waterConsumption_df.assign(InterpolatePoly5=waterConsumption_df.target.interpolate(method='polynomial', order=5))
waterConsumption_df = waterConsumption_df.assign(InterpolatePoly7=waterConsumption_df.target.interpolate(method='polynomial', order=7))
waterConsumption_df = waterConsumption_df.assign(InterpolateSpline3=waterConsumption_df.target.interpolate(method='spline', order=3))
waterConsumption_df = waterConsumption_df.assign(InterpolateSpline4=waterConsumption_df.target.interpolate(method='spline', order=4))
waterConsumption_df = waterConsumption_df.assign(InterpolateSpline5=waterConsumption_df.target.interpolate(method='spline', order=5))
# 결과에 점수를 매겨보고, 어떤 것이 더 나은지 눈으로 확인하자.
# 서로 다른 방법에 대해서 채점할 수 있는 모듈을 불러온다.
from sklearn.metrics import r2_score
results = [(method, r2_score(waterConsumption_df.reference, waterConsumption_df[method])) for method in list(waterConsumption_df)[3:]]
results_df = pd.DataFrame(np.array(results), columns=['Method', 'R_squared'])
results_df.sort_values(by='R_squared', ascending=False)
| Method | R_squared | |
|---|---|---|
| 9 | InterpolateAkima | 0.981684100149588 |
| 5 | InterpolateTime | 0.9815664478940275 |
| 8 | InterpolateSLinear | 0.9815664478940275 |
| 4 | InterpolateLinear | 0.9813215759943529 |
| 6 | InterpolateQuadratic | 0.9663474396797 |
| 12 | InterpolateSpline3 | 0.9633836918698976 |
| 7 | InterpolateCubic | 0.9633218181089737 |
| 10 | InterpolatePoly5 | 0.954157955951024 |
| 14 | InterpolateSpline5 | 0.951671359314308 |
| 11 | InterpolatePoly7 | 0.9504371542313383 |
| 13 | InterpolateSpline4 | 0.928463604189156 |
| 0 | FillMean | 0.7859894121335577 |
| 2 | RollingMean | 0.7457974578754563 |
| 1 | FillMedian | 0.7347827483233148 |
| 3 | RollingMedian | 0.6888988883243206 |
# Impute 이후의 데이터, InterpolateTime의 결과를 사용해본다. 난 InterpolateAkima가 더 잘나와서 이것을 사용한다.
final_df= waterConsumption_df[['reference', 'target', 'missing', 'InterpolateAkima' ]]
final_df.plot(style=['b-.', 'ko', 'r.', 'rx-'], figsize=(20,10));
plt.ylabel('Temperature');
plt.legend(loc='upper center', bbox_to_anchor=(0.5, 1.05),
fancybox=True, shadow=True, ncol=5, prop={'size': 14} );
한계
시계열 예측 모델은 예측, 예측에 대한 신뢰 구간 정보를 모두 만들어 낼 수 있다.
실제 관측값에 대한 위, 아래 신뢰구간
이는 예측 대비 실제 가능한 결과의 범위를 평가하거나, 모델의 동작 방법을 이해하는데 유용하게 사용된다. 예를 들어 python을 이용한 ARIMIA 모델은 모델 결과 객체를 반환하는데, forecast() 함수를 이용하면 3가지 결과 값을 보여준다.
예측에 대한 오차는 예측과 실제 데이터 간의 차이를 의미한다. 정확도를 측정하는 대표적인 지표는 RMSE, MAPE 이다.
시계열 모델은 반드시 각 값들에 대한 시간 컬럼( Unique, 입력, 적어도 예측 가능한... )을 갖고 있어야 한다. 시계열 데이터는 보통 Cleaning, Scailing, Transformation을 필요로 한다.
Frequency( 주기 ): 데이터는 모델입장에서 너무 많은 주기를 갖고 있거나, Resampling이 필요할만큼 고르지 않은 주기를 갖는 경우가 있다.
Outliers( 아주 예외적인 값, 탈 인간계...? ): 데이터는 별도로 식별 혹은 처리되어야 하는 극도로 차이나는 값들을 갖고 있을 수 있다.
Frequency:
Resampling:
Up-sampling
얼마나 세밀한 관측값이 보간법으로 사용되는지 잘 살펴보아야 한다.
resample() 함수는 pandas Dataframe으로 사용 가능하다.
Up-sampling frequency
parser = lambda x: datetime.strptime('190'+x, '%Y-%m')
shampoo_df = pd.read_csv(
os.path.join( DATA_ROOT, 'time-series-data', 'shampoo.csv' ),
header = 0, index_col = 0, parse_dates = True, squeeze = True, date_parser = parser
)
upsampled_ts = shampoo_df.resample('D').mean()
print(upsampled_ts .head(36))
Month 1901-01-01 266.0 1901-01-02 NaN 1901-01-03 NaN 1901-01-04 NaN 1901-01-05 NaN 1901-01-06 NaN 1901-01-07 NaN 1901-01-08 NaN 1901-01-09 NaN 1901-01-10 NaN 1901-01-11 NaN 1901-01-12 NaN 1901-01-13 NaN 1901-01-14 NaN 1901-01-15 NaN 1901-01-16 NaN 1901-01-17 NaN 1901-01-18 NaN 1901-01-19 NaN 1901-01-20 NaN 1901-01-21 NaN 1901-01-22 NaN 1901-01-23 NaN 1901-01-24 NaN 1901-01-25 NaN 1901-01-26 NaN 1901-01-27 NaN 1901-01-28 NaN 1901-01-29 NaN 1901-01-30 NaN 1901-01-31 NaN 1901-02-01 145.9 1901-02-02 NaN 1901-02-03 NaN 1901-02-04 NaN 1901-02-05 NaN Freq: D, Name: Sales, dtype: float64
resample()이 NaN으로 채워진 rows를 만든 것을 봤다. 다음으로 누락된 값을 새로 만들어진 주기에 맞추어 집어 넣어야 한다.interpolate 함수를 사용하면 된다. 사용 가능한 데이터 사이를 일직선으로 긋고, 그 데이터를 사용한다고 생각하면 된다.
interpolated = upsampled_ts.interpolate(method = 'linear')
interpolated.plot()
plt.show()
다른 일반적인 보간법
interpolated1 = upsampled_ts.interpolate(method = 'spline', order = 2)
interpolated1.plot()
plt.show()
print(interpolated1.head(12))
Month 1901-01-01 266.000000 1901-01-02 258.630160 1901-01-03 251.560886 1901-01-04 244.720748 1901-01-05 238.109746 1901-01-06 231.727880 1901-01-07 225.575149 1901-01-08 219.651553 1901-01-09 213.957094 1901-01-10 208.491770 1901-01-11 203.255582 1901-01-12 198.248529 Freq: D, Name: Sales, dtype: float64
Down-sampling Frequency
resample = shampoo_df.resample('Q')
quarterly_mean_sales = resample.mean()
print(quarterly_mean_sales.head())
quarterly_mean_sales.plot()
plt.show()
Month 1901-03-31 198.333333 1901-06-30 156.033333 1901-09-30 216.366667 1901-12-31 215.100000 1902-03-31 184.633333 Freq: Q-DEC, Name: Sales, dtype: float64
이제 월별 데이터를 연도별 데이터로 만들 수 있게 되었다. Down-sample을 alias로 year-end 주기로 만들되, 각 연도별 전체 판매량을 합산한 것으로 사용하자.
resample = shampoo_df.resample('A')
yearly_mean_sales = resample.sum()
print(yearly_mean_sales.head() )
yearly_mean_sales.plot()
plt.show()
Month 1901-12-31 2357.5 1902-12-31 3153.5 1903-12-31 5742.6 Freq: A-DEC, Name: Sales, dtype: float64
Outliers 데이터는 간혹 오염되거나, 매우 튀는 값을 갖는데, 이 것을 따로 식별해주거나 처리해주어야 한다.
parser = lambda x: datetime.strptime('190'+x, '%Y-%m')
shampoo_df=pd.read_csv(
os.path.join( DATA_ROOT, 'time-series-data', 'shampoo.csv' ),
header=0, index_col=0, parse_dates=True, squeeze=True, date_parser=parser
)
X = shampoo_df.values
diff = list()
for i in range(1, len(X)):
value = X[i] - X[i - 1]
diff.append(value)
plt.plot(diff)
plt.show()
from sklearn.linear_model import LinearRegression
# 선형 모델에 맞춘다.
X = [i for i in range(0, len(shampoo_df))]
X = np.reshape(X, (len(X), 1))
y = shampoo_df.values
model = LinearRegression()
model.fit(X, y)
# 추세를 계산한다.
trend = model.predict(X)
# 추세를 보자.
plt.plot(y)
plt.plot(trend)
plt.show()
# detrend
detrended = [y[i]-trend[i] for i in range(0, len(shampoo_df))]
# detrended 보기
plt.plot(detrended)
plt.show()
파란색 실제 데이터위로 오랜지색 추세선이 그려진다.
Seasonal adjustment with differencing
최저 기온 데이터를 가지고 seasonality differencing을 진행해보자.
# 월별 데이터를 Differencing을 이용하여 deseasonalize
min_temperature = pd.read_csv(
os.path.join( DATA_ROOT, 'time-series-data', 'daily-min-temperatures.csv' ),
header=0, index_col=0, parse_dates=True, squeeze=True
)
resample = min_temperature.resample('M')
monthly_mean = resample.mean()
X = min_temperature.values
diff = list()
months_in_year = 12
for i in range(months_in_year, len(monthly_mean)):
value = monthly_mean[i] - monthly_mean[i - months_in_year]
diff.append(value)
plt.plot(diff)
plt.show()
최적의 모델을 찾기 전에 보통 Moving avereage, Exponential smoothing을 사용해 본다. 모델 선택은 예측 정확도 측정방법에 의해서 결정되는데, 이를 테면
n개의 관측값이 있고, Yn은 n 시점에서의 값 Y를 의미하고, Fn은 예측된 값을 의미하고, RMSE와 MAPE는 예측 정확도를 측정할 때 사용되는 가장 대표적인 지표이다.
일별 여성 출생 데이터를 사용해보자. 데이터는 1959년 캘리포니아의 일별 여성 출생자 수를 가지고 있다. 3일 Window로 Moving average를 만들고, RMSE, MAPE로 얼마나 잘 맞는지 확인해봐라.
# MAE, MAPE 함수 선언
from sklearn.metrics import mean_squared_error
def MAE(y,yhat):
diff = np.abs(np.array(y)-np.array(yhat))
try:
mae = round(np.mean(np.fabs(diff)),3)
except:
print("Error while calculating")
mae = np.nan
return mae
def MAPE(y, yhat):
y, yhat = np.array(y), np.array(yhat)
try:
mape = round(np.mean(np.abs((y - yhat) / y)) * 100,2)
except:
print("Observed values are empty")
mape = np.nan
return mape
female_birth_series = pd.read_csv(
os.path.join( DATA_ROOT, 'time-series-data', 'daily-total-female-births.csv' ),
header=0, index_col=0, parse_dates=True, squeeze=True
)
# tail rolling average transform
rolling = female_birth_series.rolling(window = 3) # arbitrarily chosen
rolling_mean = rolling.mean()
female_birth_series.plot()
rolling_mean.plot(color = 'red')
plt.show()
# 원본 데이터와 변형된 데이터를 근접해서 살펴보자.
female_birth_series[:100].plot()
rolling_mean[:100].plot(color = 'red')
plt.show()
import numpy as np
y_df = pd.DataFrame( {'Observed':female_birth_series.values, 'Predicted':rolling_mean})
y_df .dropna(axis = 0, inplace = True)
print(y_df.tail())
rmse = np.sqrt(mean_squared_error(y_df.Observed, y_df.Predicted))
print("\n\n Accuracy measures ")
print('RMSE: %.3f' % rmse)
n = y_df.shape[0]
mae = MAE(y_df.Observed, y_df.Predicted)
print('MAE: %d' % np.float64(mae))
mape = MAPE(y_df.Observed, y_df.Predicted)
print('MAPE: %.3f' % np.float64(mape))
Observed Predicted Date 1959-12-27 37 38.333333 1959-12-28 52 41.000000 1959-12-29 48 45.666667 1959-12-30 55 51.666667 1959-12-31 50 51.000000 Accuracy measures RMSE: 5.435 MAE: 4 MAPE: 10.630
Window 크기 2, 4로 해서 다시 수행하고, 정확도가 증가되는지 아닌지 확인해보라. 샴페인 데이터를 이용해서 Seasonality 데이터를 확인해보고, 다양한 Window 크기로 Decomposing을 수행한 다음 이동 평균을 구해보아라. 그리고 RMSE, MAPE를 측정해 보아라.
1960년 Winters에 의하면 time( t + 1 )에서의 SES를 통한 예측 값은,
Parameter $\alpha$는 smoothing constant이고, 그 값은 0 ~ 1 사이에 있다. 모델이 Smoothing 상수를 하나만 사용하기 때문에 Single Exponential Smoothing이라고 부른다.
2001.01 ~ 2013.09 기간의 석유데이터를 보자.
import statsmodels.tsa.holtwinters as ets
import statsmodels.tools.eval_measures as fa
from sklearn.metrics import mean_squared_error
from statsmodels.tsa.holtwinters import SimpleExpSmoothing
def MAPE(y, yhat):
y, yhat = np.array(y), np.array(yhat)
try:
mape = round(np.sum(np.abs(yhat - y)) / np.sum(y) * 100,2)
except:
print("Observed values are empty")
mape = np.nan
return mape
petrol_df = pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'Petrol.csv' ) )
date_rng = pd.date_range(start='1/1/2001', end='30/9/2013', freq='Q')
print(date_rng)
petrol_df['TimeIndex'] = pd.DataFrame(date_rng, columns=['Quarter'])
print(petrol_df.head(3).T)
plt.plot(petrol_df.TimeIndex, petrol_df.Consumption)
plt.title('Original data before split')
plt.show()
DatetimeIndex(['2001-03-31', '2001-06-30', '2001-09-30', '2001-12-31',
'2002-03-31', '2002-06-30', '2002-09-30', '2002-12-31',
'2003-03-31', '2003-06-30', '2003-09-30', '2003-12-31',
'2004-03-31', '2004-06-30', '2004-09-30', '2004-12-31',
'2005-03-31', '2005-06-30', '2005-09-30', '2005-12-31',
'2006-03-31', '2006-06-30', '2006-09-30', '2006-12-31',
'2007-03-31', '2007-06-30', '2007-09-30', '2007-12-31',
'2008-03-31', '2008-06-30', '2008-09-30', '2008-12-31',
'2009-03-31', '2009-06-30', '2009-09-30', '2009-12-31',
'2010-03-31', '2010-06-30', '2010-09-30', '2010-12-31',
'2011-03-31', '2011-06-30', '2011-09-30', '2011-12-31',
'2012-03-31', '2012-06-30', '2012-09-30', '2012-12-31',
'2013-03-31', '2013-06-30', '2013-09-30'],
dtype='datetime64[ns]', freq='Q-DEC')
0 1 2
Year 2001 2001 2001
Quarter Q1 Q2 Q3
Consumption 14.978 11.099 10.057
TimeIndex 2001-03-31 00:00:00 2001-06-30 00:00:00 2001-09-30 00:00:00
# 학습, 테스트 데이터 생성
train = petrol_df[0:int(len(petrol_df)*0.7)]
test= petrol_df[int(len(petrol_df)*0.7):]
print("\n Training data start at \n")
print (train[train.TimeIndex == train.TimeIndex.min()],['Year','Quarter'])
print("\n Training data ends at \n")
print (train[train.TimeIndex == train.TimeIndex.max()],['Year','Quarter'])
print("\n Test data start at \n")
print (test[test.TimeIndex == test.TimeIndex.min()],['Year','Quarter'])
print("\n Test data ends at \n")
print (test[test.TimeIndex == test.TimeIndex.max()],['Year','Quarter'])
plt.plot(train.TimeIndex, train.Consumption, label = 'Train')
plt.plot(test.TimeIndex, test.Consumption, label = 'Test')
plt.legend(loc = 'best')
plt.title('Original data after split')
plt.show()
Training data start at
Year Quarter Consumption TimeIndex
0 2001 Q1 14.978 2001-03-31 ['Year', 'Quarter']
Training data ends at
Year Quarter Consumption TimeIndex
34 2009 Q3 1.98159 2009-09-30 ['Year', 'Quarter']
Test data start at
Year Quarter Consumption TimeIndex
35 2009 Q4 1.55221 2009-12-31 ['Year', 'Quarter']
Test data ends at
Year Quarter Consumption TimeIndex
50 2013 Q3 0.72823 2013-09-30 ['Year', 'Quarter']
SimpleExpSmoothing 클래스를 만들고, 학습 데이터를 입력해라.
fit() 함수는 fit smoothing_level 과 같은 fit 설정 값을 받아야 하는데, 이를 생략하거나 None 값을 넣으면, 알아서 값을 최적화한다. $\alpha$에 0.1, 0.5, 0.99와 같이 여러 값을 넣어보고, Model이 $\alpha$를 최적화 하도록 만들어야겠다.
1에 가까운 값: 빠른 학습( 가장 최근의 값들을 예측에 사용 ) 0에 가까운 값: 늦은 학습( 가장 과거의 값에 큰 영향을 받도록... )
아래 책의 Page 89을 보면된다.
# 모델 객체 생성
model = SimpleExpSmoothing(np.asarray(train['Consumption']))
# Model fitting
pred_SES = test.copy() # Have a copy of the test dataset
for alpha_value in ( 0.1, 0.5, 0.99 ):
alpha_str = "SES" + str(alpha_value)
mode_fit_i = model.fit(smoothing_level = alpha_value, optimized=False)
pred_SES[alpha_str] = mode_fit_i.forecast(len(test['Consumption']))
rmse = np.sqrt(mean_squared_error(test['Consumption'], pred_SES[alpha_str]))
mape = MAPE(test['Consumption'],pred_SES[alpha_str])
print("For alpha = %1.2f, RMSE is %3.4f MAPE is %3.2f" %(alpha_value, rmse, mape))
plt.figure(figsize=(16,8))
plt.plot(train.TimeIndex, train['Consumption'], label ='Train')
plt.plot(test.TimeIndex, test['Consumption'], label ='Test')
plt.plot(test.TimeIndex, pred_SES[alpha_str], label = alpha_str)
plt.title('Simple Exponential Smoothing with alpha ' + str(alpha_value))
plt.legend(loc='best')
plt.show()
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\holtwinters\model.py:427: FutureWarning: After 0.13 initialization must be handled at model creation warnings.warn(
For alpha = 0.10, RMSE is 3.0712 MAPE is 283.19
For alpha = 0.50, RMSE is 1.3166 MAPE is 119.40
For alpha = 0.99, RMSE is 0.9438 MAPE is 83.89
conf 값을 생략해서 최적의 $\alpha$ 값을 찾도록 하고, 그 결과를 모델에 사용할 수 있도록 해보자.
pred_opt = SimpleExpSmoothing(train['Consumption']).fit(optimized = True)
print('')
print('== Simple Exponential Smoothing ')
print('')
print('')
print('Smoothing Level', np.round(pred_opt.params['smoothing_level'], 4))
print('Initial Level', np.round(pred_opt.params['initial_level'], 4))
print('')
y_pred_opt = pred_opt.forecast(steps = 16)
df_pred_opt = pd.DataFrame({'Y_hat':y_pred_opt,'Y':test['Consumption'].values})
rmse_opt = np.sqrt(mean_squared_error(test['Consumption'], y_pred_opt))
mape_opt = MAPE(test['Consumption'], y_pred_opt)
alpha_value = np.round(pred_opt.params['smoothing_level'], 4)
print("For alpha = %1.2f, RMSE is %3.4f MAPE is %3.2f" %(alpha_value, rmse_opt, mape_opt))
plt.figure(figsize=(16,8))
plt.plot(train.TimeIndex, train['Consumption'], label = 'Train')
plt.plot(test.TimeIndex, test['Consumption'], label = 'Test')
plt.plot(test.TimeIndex, y_pred_opt, label = 'SES_OPT')
plt.title('Simple Exponential Smoothing with alpha ' + str(alpha_value))
plt.legend(loc='best')
plt.show()
print(df_pred_opt.head().T)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\holtwinters\model.py:427: FutureWarning: After 0.13 initialization must be handled at model creation warnings.warn(
== Simple Exponential Smoothing Smoothing Level 1.0 Initial Level 14.978 For alpha = 1.00, RMSE is 0.9386 MAPE is 83.39
35 36 37 38 39 Y_hat 1.98159 1.98159 1.98159 1.98159 1.98159 Y 1.55221 1.57034 1.08986 1.05821 1.39665
최적의 $\alpha$ 값을 확인하고, RMSW, MAPE도 다른 $\alpha$ 값보다 작은 것을 확인했다
다른 방법도 해 볼 수 있다.
Double Exponential Smoothing는 예측을 위해서 2개 수식을 사용한다. 하나는 단기 예측 평균 값 혹은 Level이고, 나머지 하나는 추세 탐지를 위한 것이다.
Intercept or Level equation, $L_t$ is given by: $L_t = {\alpha}{Y_t} + (1 - \alpha)F_t$
Trend equation is given by $T_t = {\beta}{(L_t - L_{t-1})} + (1 - \beta)T_{t-1}$
Here, $\alpha$ and $\beta$ 들은 각각 Level, Trend를 위한 Smoothing 상수이다.
t + 1 시점에서의 예측은 아래와 같다.
아래와 같은 경우를 생각해보자.
from statsmodels.tsa.holtwinters import Holt
model = Holt(np.asarray(train['Consumption']))
model_fit = model.fit()
print('')
print('==Holt model Exponential Smoothing Parameters ==')
print('')
alpha_value = np.round(model_fit.params['smoothing_level'], 4)
print('Smoothing Level', alpha_value )
print('Smoothing Trend', np.round(model_fit.params['smoothing_trend'], 4))
print('Initial Level', np.round(model_fit.params['initial_level'], 4))
print('')
==Holt model Exponential Smoothing Parameters == Smoothing Level 1.0 Smoothing Trend 0.0 Initial Level 15.3602
Pred_Holt = test.copy()
Pred_Holt['Opt'] = model_fit.forecast(len(test['Consumption']))
plt.figure(figsize=(16,8))
plt.plot(train['Consumption'], label='Train')
plt.plot(test['Consumption'], label='Test')
plt.plot(Pred_Holt['Opt'], label='HoltOpt')
plt.legend(loc='best')
plt.show()
df_pred_opt = pd.DataFrame({'Y_hat':Pred_Holt['Opt'] ,'Y':test['Consumption'].values})
rmse_opt = np.sqrt(mean_squared_error(df_pred_opt.Y, df_pred_opt.Y_hat))
mape_opt = MAPE(df_pred_opt.Y, df_pred_opt.Y_hat)
print("For alpha = %1.2f, RMSE is %3.4f MAPE is %3.2f" %(alpha_value, rmse_opt, mape_opt))
For alpha = 1.00, RMSE is 2.8192 MAPE is 217.85
# 모델 파라미터는?
print(model_fit.params)
{'smoothing_level': 0.9999999850988388, 'smoothing_trend': 0.0, 'smoothing_seasonal': nan, 'damping_trend': nan, 'initial_level': 15.360248414159342, 'initial_trend': -0.382239607547084, 'initial_seasons': array([], dtype=float64), 'use_boxcox': False, 'lamda': None, 'remove_bias': False}
동일한 석유 데이터 사용했을 때, SES, HOLT 모델의 결과는 다르다. ( 원래 자료에서는 Holt가 더 좋음. )
| Model | RMSE | MAPE |
|---|---|---|
| SES | 0.9386 | 83.39 |
| Holt | 2.8192 | 217.85 |
아래 방법으로도 시도해 볼 수 있다.
optimized = False
동일한 모델링을 Training data로 해봐라.
from statsmodels.tsa.holtwinters import ExponentialSmoothing
pred1 = ExponentialSmoothing(
np.asarray(train['Consumption']),
trend='additive', damped=False, seasonal='additive',
seasonal_periods = 12
).fit() #[:'2017-01-01']
print('')
print('== Holt-Winters Additive ETS(A,A,A) Parameters ==')
print('')
alpha_value = np.round(pred1.params['smoothing_level'], 4)
print('Smoothing Level: ', alpha_value)
print('Smoothing Trend: ', np.round(pred1.params['smoothing_trend'], 4))
print('Smoothing Seasonal: ', np.round(pred1.params['smoothing_seasonal'], 4))
print('Initial Level: ', np.round(pred1.params['initial_level'], 4))
print('Initial Trend: ', np.round(pred1.params['initial_trend'], 4))
print('Initial Seasons: ', np.round(pred1.params['initial_seasons'], 4))
print('')
### Forecast for next 16 months
y_pred1 = pred1.forecast(steps = 16)
df_pred1 = pd.DataFrame({'Y_hat':y_pred1,'Y':test['Consumption']})
print(df_pred1)
<ipython-input-117-8ced8d035bf7>:3: FutureWarning: the 'damped'' keyword is deprecated, use 'damped_trend' instead pred1 = ExponentialSmoothing( C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\holtwinters\model.py:427: FutureWarning: After 0.13 initialization must be handled at model creation warnings.warn(
== Holt-Winters Additive ETS(A,A,A) Parameters ==
Smoothing Level: 1.0
Smoothing Trend: 0.0
Smoothing Seasonal: 0.0
Initial Level: 12.8762
Initial Trend: -0.3689
Initial Seasons: [ 2.4708 0.0845 -0.3991 0.169 0.6992 -0.355 0.4698 1.9918 2.7108
1.7961 2.0184 2.5988]
Y_hat Y
35 2.193087 1.55221
36 1.696143 1.57034
37 -1.059051 1.08986
38 -1.911502 1.05821
39 -1.712360 1.39665
40 -1.551110 1.32161
41 -2.974147 0.96034
42 -2.518340 0.87553
43 -1.365209 1.07530
44 -1.015152 1.30285
45 -2.298785 0.88939
46 -2.445392 0.88818
47 -2.233894 0.99804
48 -2.730839 0.84120
49 -5.486033 0.74032
50 -6.338483 0.72823
### Plot
fig2, ax = plt.subplots()
ax.plot(df_pred1.Y, label='Original')
ax.plot(df_pred1.Y_hat, label='Predicted')
plt.legend(loc='upper left')
plt.title('Holt-Winters Additive ETS(A,A,A) Method 1')
plt.ylabel('Qty')
plt.xlabel('Date')
plt.show()
# 위 그래프가 원본과 완전히 벗어남. HOLT 못 써먹을 듯...인가??? 아니면 본 데이터에 맞는 모델을 저자가 썼었나...?
rmse = np.sqrt(mean_squared_error(df_pred1.Y, df_pred1.Y_hat))
mape = MAPE(df_pred1.Y, df_pred1.Y_hat)
print("For alpha = %1.2f, RMSE is %3.4f MAPE is %3.2f" %(alpha_value, rmse, mape))
For alpha = 1.00, RMSE is 3.5646 MAPE is 292.53
print(pred1.params)
{'smoothing_level': 0.9999999850988388, 'smoothing_trend': 0.0, 'smoothing_seasonal': 8.597195499160461e-10, 'damping_trend': nan, 'initial_level': 12.876160591070663, 'initial_trend': -0.36891513981272395, 'initial_seasons': array([ 2.47075063, 0.084472 , -0.3990635 , 0.16899335, 0.69915859,
-0.35496347, 0.46975926, 1.99180547, 2.71077743, 1.79605932,
2.01836774, 2.59878007]), 'use_boxcox': False, 'lamda': None, 'remove_bias': False}
Holt-Winters Additive ETS( A, A, A ) 모델의 MAPE: 128.21( 원본 ) 292.53( 실제 ) Holt-Winters Additive ETS( A, A, N ) 모델의 MAPE: 63.01( 원본 ) 217.85( 실제 )
실제로 ETS( A, A, A )의 MAPE가 높기는 한데... 영 아닌것 같다.... 너무 안 좋다...
항공 여객 데이터로 해보자.
AirPax = pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'AirPax.csv' ) )
date_rng = pd.date_range(start='1/1/1949', end='31/12/1960', freq='M')
print(date_rng)
AirPax['TimeIndex'] = pd.DataFrame(date_rng, columns=['Month'])
print(AirPax.head())
DatetimeIndex(['1949-01-31', '1949-02-28', '1949-03-31', '1949-04-30',
'1949-05-31', '1949-06-30', '1949-07-31', '1949-08-31',
'1949-09-30', '1949-10-31',
...
'1960-03-31', '1960-04-30', '1960-05-31', '1960-06-30',
'1960-07-31', '1960-08-31', '1960-09-30', '1960-10-31',
'1960-11-30', '1960-12-31'],
dtype='datetime64[ns]', length=144, freq='M')
Year Month Passenger TimeIndex
0 1949 Jan 112 1949-01-31
1 1949 Feb 118 1949-02-28
2 1949 Mar 132 1949-03-31
3 1949 Apr 129 1949-04-30
4 1949 May 121 1949-05-31
# 학습, 테스트 데이터 셋 만들기. 7:3으로 자름.
train = AirPax[0:int(len(AirPax)*0.7)]
test= AirPax[int(len(AirPax)*0.7):]
print("\n Training data start at \n")
print (train[train.TimeIndex == train.TimeIndex.min()],['Year','Month'])
print("\n Training data ends at \n")
print (train[train.TimeIndex == train.TimeIndex.max()],['Year','Month'])
print("\n Test data start at \n")
print (test[test.TimeIndex == test.TimeIndex.min()],['Year','Month'])
print("\n Test data ends at \n")
print (test[test.TimeIndex == test.TimeIndex.max()],['Year','Month'])
plt.plot(train.TimeIndex, train.Passenger, label = 'Train')
plt.plot(test.TimeIndex, test.Passenger, label = 'Test')
plt.legend(loc = 'best')
plt.title('Original data after split')
plt.show()
Training data start at
Year Month Passenger TimeIndex
0 1949 Jan 112 1949-01-31 ['Year', 'Month']
Training data ends at
Year Month Passenger TimeIndex
99 1957 Apr 348 1957-04-30 ['Year', 'Month']
Test data start at
Year Month Passenger TimeIndex
100 1957 May 355 1957-05-31 ['Year', 'Month']
Test data ends at
Year Month Passenger TimeIndex
143 1960 Dec 432 1960-12-31 ['Year', 'Month']
pred = ExponentialSmoothing(np.asarray(train['Passenger']),
seasonal_periods=12 ,seasonal='add').fit(optimized=True)
print(pred.params)
print('')
print('== Holt-Winters Additive ETS(A,A,A) Parameters ==')
print('')
alpha_value = np.round(pred.params['smoothing_level'], 4)
print('Smoothing Level: ', alpha_value)
print('Smoothing Trend: ', np.round(pred.params['smoothing_trend'], 4))
print('Smoothing Seasonal: ', np.round(pred.params['smoothing_seasonal'], 4))
print('Initial Level: ', np.round(pred.params['initial_level'], 4))
print('Initial Trend: ', np.round(pred.params['initial_trend'], 4))
print('Initial Seasons: ', np.round(pred.params['initial_seasons'], 4))
print('')
{'smoothing_level': 0.9999988342354306, 'smoothing_trend': nan, 'smoothing_seasonal': 1.0606589788585584e-06, 'damping_trend': nan, 'initial_level': 202.28028252920663, 'initial_trend': nan, 'initial_seasons': array([ -90.2711059 , -93.97460735, -64.20906334, -72.33395774,
-74.60152474, -48.87540266, -24.39321575, -28.41077563,
-58.94271216, -89.33333766, -116.35329845, -93.00101449]), 'use_boxcox': False, 'lamda': None, 'remove_bias': False}
== Holt-Winters Additive ETS(A,A,A) Parameters ==
Smoothing Level: 1.0
Smoothing Trend: nan
Smoothing Seasonal: 0.0
Initial Level: 202.2803
Initial Trend: nan
Initial Seasons: [ -90.2711 -93.9746 -64.2091 -72.334 -74.6015 -48.8754 -24.3932
-28.4108 -58.9427 -89.3333 -116.3533 -93.001 ]
pred_HoltW = test.copy()
pred_HoltW['HoltW'] = model_fit.forecast(len(test['Passenger']))
plt.figure(figsize=(16,8))
plt.plot(train['Passenger'], label='Train')
plt.plot(test['Passenger'], label='Test')
plt.plot(pred_HoltW['HoltW'], label='HoltWinters')
plt.title('Holt-Winters Additive ETS(A,A,A) Parameters:\n alpha = ' +
str(alpha_value) + ' Beta:' +
str(np.round(pred.params['smoothing_trend'], 4)) +
' Gamma: ' + str(np.round(pred.params['smoothing_seasonal'], 4)))
plt.legend(loc='best')
plt.show()
df_pred_opt = pd.DataFrame({'Y_hat':pred_HoltW['HoltW'] ,'Y':test['Passenger'].values})
rmse_opt = np.sqrt(mean_squared_error(df_pred_opt.Y, df_pred_opt.Y_hat))
mape_opt = MAPE(df_pred_opt.Y, df_pred_opt.Y_hat)
print("For alpha = %1.2f, RMSE is %3.4f MAPE is %3.2f" %(alpha_value, rmse_opt, mape_opt))
# 쓰지 말아야 겠다... 이건 뭐... 왜 이렇게 나오는지도 모르겠다...
For alpha = 1.00, RMSE is 434.9059 MAPE is 101.57
항공 여객 데이터, ETS(A, A, M) model을 이용해서 지난 12개월을 예측해보자.
pred = ExponentialSmoothing(
np.asarray(train['Passenger']),
seasonal_periods=12 ,seasonal='multiplicative'
).fit(optimized=True)
print('')
print('== Holt-Winters Additive ETS(A,A,A) Parameters ==')
print('')
alpha_value = np.round(pred.params['smoothing_level'], 4)
print('Smoothing Level: ', alpha_value)
print('Smoothing Trend: ', np.round(pred.params['smoothing_trend'], 4))
print('Smoothing Seasonal: ', np.round(pred.params['smoothing_seasonal'], 4))
print('Initial Level: ', np.round(pred.params['initial_level'], 4))
print('Initial Trend: ', np.round(pred.params['initial_trend'], 4))
print('Initial Seasons: ', np.round(pred.params['initial_seasons'], 4))
print('')
== Holt-Winters Additive ETS(A,A,A) Parameters == Smoothing Level: 0.8806 Smoothing Trend: nan Smoothing Seasonal: 0.0 Initial Level: 141.0948 Initial Trend: nan Initial Seasons: [0.801 0.7769 0.8975 0.8694 0.8627 0.9751 1.0809 1.0605 0.9327 0.8087 0.6991 0.7893]
pred_HoltW = test.copy()
pred_HoltW['HoltWM'] = pred.forecast(len(test['Passenger']))
plt.figure(figsize=(16,8))
plt.plot(train['Passenger'], label='Train')
plt.plot(test['Passenger'], label='Test')
plt.plot(pred_HoltW['HoltWM'], label='HoltWinters')
plt.title('Holt-Winters Multiplicative ETS(A,A,M) Parameters:\n alpha = ' +
str(alpha_value) + ' Beta:' +
str(np.round(pred.params['smoothing_trend'], 4)) +
' Gamma: ' + str(np.round(pred.params['smoothing_seasonal'], 4)))
plt.legend(loc='best')
plt.show()
# 이건 좀 비슷하게라도 가네....
df_pred_opt = pd.DataFrame({'Y_hat':pred_HoltW['HoltWM'] ,'Y':test['Passenger'].values})
rmse_opt = np.sqrt(mean_squared_error(df_pred_opt.Y, df_pred_opt.Y_hat))
mape_opt = MAPE(df_pred_opt.Y, df_pred_opt.Y_hat)
print("For alpha = %1.2f, RMSE is %3.4f MAPE is %3.2f" %(alpha_value, rmse_opt, mape_opt))
For alpha = 0.88, RMSE is 83.4364 MAPE is 16.14
MAPE of ETS(A,A,M) model is 16.14 and it is lesser than 101.57, MAPE of ETS(A,A,A) model 라고는 하는데... ETS( A, A, A )는 박살난 결과였기에 어디 비교할 대상도 안됨.
반면, Holt-Winter - Additive Trend and Multiplicative Seasonality model이 Holt-Winter - Additive model 보다 낫다.
샴푸 데이터 써보자. 3년동안의 월간 샴푸 판매량이 나와 있다.총 36개의 관측값이 있다. The original dataset is credited to Makridakis, Wheelwright and Hyndman (1998). 전체 데이터는 아래 URL에 있다.
def parser(x):
return pd.datetime.strptime('190'+x, '%Y-%m')
shampoo_df = pd.read_csv(
os.path.join( DATA_ROOT, 'time-series-data', 'shampoo.csv' ),
header=0, parse_dates=True, squeeze=True, date_parser=parser
)
print(shampoo_df.head())
shampoo_df.plot()
Month Sales 0 1-01 266.0 1 1-02 145.9 2 1-03 183.1 3 1-04 119.3 4 1-05 180.3
<AxesSubplot:>
# 학습, 테스트 데이터 분리
train = shampoo_df[0:int(len(shampoo_df)*0.7)]
test = shampoo_df[int(len(shampoo_df)*0.7):]
train['Sales'].plot(figsize=(15,8), title= 'Monthly Sales', fontsize=14)
test['Sales'].plot(figsize=(15,8), title= 'Monthly Sales', fontsize=14)
shampoo_df.head()
| Month | Sales | |
|---|---|---|
| 0 | 1-01 | 266.0 |
| 1 | 1-02 | 145.9 |
| 2 | 1-03 | 183.1 |
| 3 | 1-04 | 119.3 |
| 4 | 1-05 | 180.3 |
# 1. 시간 회귀
shampoo_df1 = shampoo_df.copy() # 복사본 하나 만들고
time = [i+1 for i in range(len(shampoo_df))]
shampoo_df1['time'] = time
monthDf = shampoo_df1[['Month']]
shampoo_df1.drop('Month', axis=1, inplace=True)
shampoo_df1.head(2)
| Sales | time | |
|---|---|---|
| 0 | 266.0 | 1 |
| 1 | 145.9 | 2 |
train=shampoo_df1[0:int(len(shampoo_df1)*0.7)]
test=shampoo_df1[int(len(shampoo_df1)*0.7):]
x_train = train.drop('Sales', axis=1)
x_test = test.drop('Sales', axis=1)
y_train = train[['Sales']]
y_test = test[['Sales']]
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(x_train, y_train)
LinearRegression()
predictions = model.predict(x_test)
y_test['RegOnTime'] = predictions
plt.figure(figsize=(16,8))
plt.plot( train['Sales'], label='Train')
plt.plot(test['Sales'], label='Test')
plt.plot(y_test['RegOnTime'], label='Regression On Time')
plt.legend(loc='best')
<matplotlib.legend.Legend at 0x2126a0b5b80>
from math import sqrt
rmse = sqrt(mean_squared_error(test.Sales, y_test.RegOnTime))
rmse = round(rmse, 3)
mape = MAPE(test.Sales, y_test.RegOnTime)
print("For RegressionOnTime, RMSE is %3.3f MAPE is %3.2f" %(rmse, mape))
For RegressionOnTime, RMSE is 164.563 MAPE is 27.94
resultsDf = pd.DataFrame({'Method':['RegressionOnTime'], 'rmse': [rmse], 'mape' : [mape]})
resultsDf
| Method | rmse | mape | |
|---|---|---|---|
| 0 | RegressionOnTime | 164.563 | 27.94 |
time = [i+1 for i in range(len(shampoo_df))]
shampoo_df1 = shampoo_df.copy()
shampoo_df1['time'] = time
print(shampoo_df1.head())
print(shampoo_df1.shape[0])
monthSeasonality = ['m1', 'm2', 'm3', 'm4', 'm5', 'm6', 'm7', 'm8', 'm9', 'm10', 'm11', 'm12']
Month Sales time 0 1-01 266.0 1 1 1-02 145.9 2 2 1-03 183.1 3 3 1-04 119.3 4 4 1-05 180.3 5 36
shampoo_df1['monthSeasonality'] = monthSeasonality * 3
shampoo_df1.head()
| Month | Sales | time | monthSeasonality | |
|---|---|---|---|---|
| 0 | 1-01 | 266.0 | 1 | m1 |
| 1 | 1-02 | 145.9 | 2 | m2 |
| 2 | 1-03 | 183.1 | 3 | m3 |
| 3 | 1-04 | 119.3 | 4 | m4 |
| 4 | 1-05 | 180.3 | 5 | m5 |
monthDf = shampoo_df1[['Month']]
shampoo_df1.drop('Month', axis=1, inplace=True)
shampoo_df1Complete = pd.get_dummies(shampoo_df1, drop_first=True)
shampoo_df1Complete.head(2).T
| 0 | 1 | |
|---|---|---|
| Sales | 266.0 | 145.9 |
| time | 1.0 | 2.0 |
| monthSeasonality_m10 | 0.0 | 0.0 |
| monthSeasonality_m11 | 0.0 | 0.0 |
| monthSeasonality_m12 | 0.0 | 0.0 |
| monthSeasonality_m2 | 0.0 | 1.0 |
| monthSeasonality_m3 | 0.0 | 0.0 |
| monthSeasonality_m4 | 0.0 | 0.0 |
| monthSeasonality_m5 | 0.0 | 0.0 |
| monthSeasonality_m6 | 0.0 | 0.0 |
| monthSeasonality_m7 | 0.0 | 0.0 |
| monthSeasonality_m8 | 0.0 | 0.0 |
| monthSeasonality_m9 | 0.0 | 0.0 |
#Creating train and test set
train=shampoo_df1Complete[0:int(len(shampoo_df1Complete)*0.7)]
test=shampoo_df1Complete[int(len(shampoo_df1Complete)*0.7):]
x_train = train.drop('Sales', axis=1)
x_test = test.drop('Sales', axis=1)
y_train = train[['Sales']]
y_test = test[['Sales']]
model = LinearRegression()
model.fit(x_train, y_train)
predictions = model.predict(x_test)
y_test['RegOnTimeSeasonal'] = predictions
<ipython-input-159-e296f75f19aa>:4: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy y_test['RegOnTimeSeasonal'] = predictions
plt.figure(figsize=(16,8))
plt.plot( train['Sales'], label='Train')
plt.plot(test['Sales'], label='Test')
plt.plot(y_test['RegOnTimeSeasonal'], label='Regression On Time With Seasonal Components')
plt.legend(loc='best')
<matplotlib.legend.Legend at 0x21269f91f10>
rmse = sqrt(mean_squared_error(test.Sales, y_test.RegOnTimeSeasonal))
rmse = round(rmse, 3)
mape = MAPE(test.Sales, y_test.RegOnTimeSeasonal)
print("For RegOnTimeSeasonal, RMSE is %3.3f MAPE is %3.2f" %(rmse, mape))
For RegOnTimeSeasonal, RMSE is 185.864 MAPE is 32.59
tempResultsDf = pd.DataFrame({'Method':['RegressionOnTimeSeasonal'], 'rmse': [rmse], 'mape' : [mape]})
resultsDf = pd.concat([resultsDf, tempResultsDf])
resultsDf
| Method | rmse | mape | |
|---|---|---|---|
| 0 | RegressionOnTime | 164.563 | 27.94 |
| 0 | RegressionOnTimeSeasonal | 185.864 | 32.59 |
# 이걸 모델이라고 해도 되나...
dd= np.asarray(train.Sales)
y_hat = test.copy()
y_hat['naive'] = dd[len(dd)-1]
plt.figure(figsize=(12,8))
plt.plot(train.index, train['Sales'], label='Train')
plt.plot(test.index,test['Sales'], label='Test')
plt.plot(y_hat.index,y_hat['naive'], label='Naive Forecast')
plt.legend(loc='best')
plt.title("Naive Forecast")
Text(0.5, 1.0, 'Naive Forecast')
rmse = sqrt(mean_squared_error(test.Sales, y_hat.naive))
rmse = round(rmse, 3)
mape = MAPE(test.Sales, y_hat.naive)
print("For Naive model, RMSE is %3.3f MAPE is %3.2f" %(rmse, mape))
For Naive model, RMSE is 186.469 MAPE is 31.72
tempResultsDf = pd.DataFrame({'Method':['Naive model'], 'rmse': [rmse], 'mape' : [mape]})
resultsDf = pd.concat([resultsDf, tempResultsDf])
resultsDf
| Method | rmse | mape | |
|---|---|---|---|
| 0 | RegressionOnTime | 164.563 | 27.94 |
| 0 | RegressionOnTimeSeasonal | 185.864 | 32.59 |
| 0 | Naive model | 186.469 | 31.72 |
RMSE, MAPE 값으로 보았을 때... Naive 방법, Regression with seasonal components 방법은 데이터가 높은 변동성을 가진 경우 사용하지 않는게 좋아보인다. Naive 방법은 안정적인 데이터셋인 경우 유효하다.
다른 방버으로 성능을 올릴 수 있을 것 같다. 다른 테크닉을 알아보자.( 기대 안되기 시작함... )
y_hat_avg = test.copy()
y_hat_avg['avg_forecast'] = train['Sales'].mean()
plt.figure(figsize=(12,8))
plt.plot(train['Sales'], label='Train')
plt.plot(test['Sales'], label='Test')
plt.plot(y_hat_avg['avg_forecast'], label='Average Forecast')
plt.legend(loc='best')
<matplotlib.legend.Legend at 0x2126a066c10>
rmse = sqrt(mean_squared_error(test.Sales, y_hat_avg.avg_forecast))
rmse = round(rmse, 3)
mape = MAPE(test.Sales, y_hat_avg.avg_forecast)
print("For Simple Average model, RMSE is %3.3f MAPE is %3.2f" %(rmse, mape))
For Simple Average model, RMSE is 279.196 MAPE is 52.35
tempResultsDf = pd.DataFrame({'Method':['Simple Average'], 'rmse': [rmse], 'mape' : [mape]})
resultsDf = pd.concat([resultsDf, tempResultsDf])
resultsDf
| Method | rmse | mape | |
|---|---|---|---|
| 0 | RegressionOnTime | 164.563 | 27.94 |
| 0 | RegressionOnTimeSeasonal | 185.864 | 32.59 |
| 0 | Naive model | 186.469 | 31.72 |
| 0 | Simple Average | 279.196 | 52.35 |
점수로 보건데 모델이 전혀 좋아지지 않았다. 반면, 각 주기의 데이터가 상수 값으로 유지 되었다면 본 점수는 좋았을 것이다.( 당연한 이야기..? ) Naive 모델의 결과 값이 더 높게 나왔다고 해서 Naive 모델이 모든 데이터 셋에서 단순 평균 모델보다 더 성능이 낫다는 것은 아니다.다른 모델이 사용해서 성능을 향상 시킬수 있는지 확인해 보자.
shampoo_df1 = shampoo_df.copy()
shampoo_df1['moving_avg_forecast_4'] = shampoo_df['Sales'].rolling(4).mean()
shampoo_df1['moving_avg_forecast_6'] = shampoo_df['Sales'].rolling(6).mean()
shampoo_df1['moving_avg_forecast_8'] = shampoo_df['Sales'].rolling(8).mean()
shampoo_df1['moving_avg_forecast_12'] = shampoo_df['Sales'].rolling(12).mean()
cols = ['moving_avg_forecast_4','moving_avg_forecast_6','moving_avg_forecast_8','moving_avg_forecast_12']
train=shampoo_df1[0:int(len(shampoo_df1)*0.7)]
test=shampoo_df1[int(len(shampoo_df1)*0.7):]
y_hat_avg = test.copy()
for col_name in cols:
plt.figure(figsize=(16,8))
plt.plot(train['Sales'], label='Train')
plt.plot(test['Sales'], label='Test')
plt.plot(y_hat_avg[col_name], label = col_name)
plt.legend(loc = 'best')
rmse = sqrt(mean_squared_error(test.Sales, y_hat_avg[col_name]))
rmse = round(rmse, 3)
mape = MAPE(test.Sales, y_hat_avg[col_name])
print("For Simple Average model, %s RMSE is %3.3f MAPE is %3.2f" %(col_name, rmse, mape))
tempResultsDf = pd.DataFrame({'Method':[col_name], 'rmse': [rmse], 'mape' : [mape]})
resultsDf = pd.concat([resultsDf, tempResultsDf])
For Simple Average model, moving_avg_forecast_4 RMSE is 75.514 MAPE is 13.05 For Simple Average model, moving_avg_forecast_6 RMSE is 88.576 MAPE is 14.62 For Simple Average model, moving_avg_forecast_8 RMSE is 105.158 MAPE is 17.03 For Simple Average model, moving_avg_forecast_12 RMSE is 134.071 MAPE is 22.79
print(resultsDf)
Method rmse mape 0 RegressionOnTime 164.563 27.94 0 RegressionOnTimeSeasonal 185.864 32.59 0 Naive model 186.469 31.72 0 Simple Average 279.196 52.35 0 moving_avg_forecast_4 75.514 13.05 0 moving_avg_forecast_6 88.576 14.62 0 moving_avg_forecast_8 105.158 17.03 0 moving_avg_forecast_12 134.071 22.79
아직 까지는 윈도우 4로 이동 평균 돌린 모델의 MAPE, RMSE 값이 가장 작다. 다른 모델도 보자.
from statsmodels.tsa.api import ExponentialSmoothing, SimpleExpSmoothing, Holt
# 모델 선언
model = SimpleExpSmoothing(train['Sales'])
model_fit = model.fit(optimized = True)
print('')
print('== Simple Exponential Smoothing ')
print('')
print('')
print('Smoothing Level', np.round(model_fit.params['smoothing_level'], 4))
print('Initial Level', np.round(model_fit.params['initial_level'], 4))
print('')
== Simple Exponential Smoothing Smoothing Level 0.2955 Initial Level 196.3512
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\holtwinters\model.py:427: FutureWarning: After 0.13 initialization must be handled at model creation warnings.warn(
y_hat_avg['SES']= model_fit.forecast(len(test['Sales']))
alpha_value = np.round(model_fit.params['smoothing_level'], 4)
plt.figure(figsize=(16,8))
plt.plot(train.index, train['Sales'], label = 'Train')
plt.plot(test.index, test['Sales'], label = 'Test')
plt.plot(test.index, y_hat_avg.SES, label = 'SES_OPT')
plt.title('Simple Exponential Smoothing with alpha ' + str(alpha_value))
plt.legend(loc='best')
plt.show()
rmse_opt= np.sqrt(mean_squared_error(test['Sales'], y_hat_avg.SES))
mape_opt= MAPE(test['Sales'], y_hat_avg.SES)
print("For alpha = %1.2f, RMSE is %3.4f MAPE is %3.2f" %(alpha_value, rmse_opt, mape_opt))
tempResultsDf = pd.DataFrame({'Method': 'SES', 'rmse': [rmse_opt], 'mape' : [mape_opt]})
resultsDf = pd.concat([resultsDf, tempResultsDf])
print(resultsDf) # 기대했던 대로... 쓰면 안된다...
For alpha = 0.30, RMSE is 203.7066 MAPE is 35.18
Method rmse mape
0 RegressionOnTime 164.563000 27.94
0 RegressionOnTimeSeasonal 185.864000 32.59
0 Naive model 186.469000 31.72
0 Simple Average 279.196000 52.35
0 moving_avg_forecast_4 75.514000 13.05
0 moving_avg_forecast_6 88.576000 14.62
0 moving_avg_forecast_8 105.158000 17.03
0 moving_avg_forecast_12 134.071000 22.79
0 SES 203.706554 35.18
0 SES 203.706554 35.18
import statsmodels.api as sm
y_hat_avg = test.copy()
model_fit = Holt(np.asarray(train['Sales'])).fit()
y_hat_avg['Holt_linear'] = model_fit.forecast(len(test))
print('')
print('==Holt model Exponential Smoothing Parameters ==')
print('')
alpha_value = np.round(model_fit.params['smoothing_level'], 4)
beta_value = np.round(model_fit.params['smoothing_trend'], 4)
print('Smoothing Level', alpha_value )
print('Smoothing Trend', beta_value)
print('Initial Level', np.round(model_fit.params['initial_level'], 4))
print('')
==Holt model Exponential Smoothing Parameters == Smoothing Level 0.0625 Smoothing Trend 0.0364 Initial Level 146.9359
plt.figure(figsize=(16,8))
plt.plot(train['Sales'], label='Train')
plt.plot(test['Sales'], label='Test')
plt.plot(y_hat_avg['Holt_linear'], label='Holt_linear')
plt.legend(loc='best')
plt.show()
rmse_opt = np.sqrt(mean_squared_error(test['Sales'], y_hat_avg['Holt_linear']))
mape_opt = MAPE(test['Sales'], y_hat_avg['Holt_linear'])
print("For alpha = %1.2f, RMSE is %3.4f MAPE is %3.2f" %(alpha_value, rmse_opt, mape_opt))
# 이쯤 되면 본 데이터에 이 모델을 왜 쓰는지 잘 이해가 안간다;;; 잘 맞는 데이터를 가져왔으면 모를까...( 심지어 원본에서도 안맞는다. )
For alpha = 0.06, RMSE is 164.8712 MAPE is 28.00
tempResultsDf = pd.DataFrame({'Method': 'Holt_linear', 'rmse': [rmse_opt], 'mape' : [mape_opt]})
resultsDf = pd.concat([resultsDf, tempResultsDf])
print(resultsDf)
Method rmse mape 0 RegressionOnTime 164.563000 27.94 0 RegressionOnTimeSeasonal 185.864000 32.59 0 Naive model 186.469000 31.72 0 Simple Average 279.196000 52.35 0 moving_avg_forecast_4 75.514000 13.05 0 moving_avg_forecast_6 88.576000 14.62 0 moving_avg_forecast_8 105.158000 17.03 0 moving_avg_forecast_12 134.071000 22.79 0 SES 203.706554 35.18 0 SES 203.706554 35.18 0 Holt_linear 164.871190 28.00
SES 보다는 조금 낫다. 이동 평균이랑은 비교하면 안된다.
y_hat_avg = test.copy()
model_fit = ExponentialSmoothing(np.asarray(train['Sales']) ,seasonal_periods = 12 ,trend='add', seasonal='add').fit()
y_hat_avg['Holt_Winter'] = model_fit.forecast(len(test))
print('')
print('== Holt-Winters Additive ETS(A,A,A) Parameters ==')
print('')
alpha_value = np.round(pred.params['smoothing_level'], 4)
beta_value = np.round(model_fit.params['smoothing_trend'], 4)
gamma_value = np.round(model_fit.params['smoothing_seasonal'], 4)
print('Smoothing Level: ', alpha_value)
print('Smoothing Trend: ', beta_value)
print('Smoothing Seasonal: ', gamma_value)
print('Initial Level: ', np.round(model_fit.params['initial_level'], 4))
print('Initial Trend: ', np.round(model_fit.params['initial_trend'], 4))
print('Initial Seasons: ', np.round(model_fit.params['initial_seasons'], 4))
print('')
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\holtwinters\model.py:427: FutureWarning: After 0.13 initialization must be handled at model creation warnings.warn(
== Holt-Winters Additive ETS(A,A,A) Parameters == Smoothing Level: 0.8806 Smoothing Trend: 0.0 Smoothing Seasonal: 0.0002 Initial Level: 234.0932 Initial Trend: 5.5282 Initial Seasons: [ -40.9379 -130.1421 -86.5721 -92.372 -108.8568 -72.036 -76.5679 -46.4252 -74.7764 -49.6152 -27.9393 -68.7575]
plt.figure(figsize=(16,8))
plt.plot( train['Sales'], label='Train')
plt.plot(test['Sales'], label='Test')
plt.plot(y_hat_avg['Holt_Winter'], label='Holt_Winter')
plt.title('Holt-Winters Multiplicative ETS(A,A,M) Parameters:\n alpha = ' +
str(alpha_value) + ' Beta:' +
str(beta_value) +
' Gamma: ' + str(gamma_value))
plt.legend(loc='best')
<matplotlib.legend.Legend at 0x2126b935910>
rmse_opt= np.sqrt(mean_squared_error(test['Sales'], y_hat_avg['Holt_Winter']))
mape_opt= MAPE(test['Sales'], y_hat_avg['Holt_Winter'])
print("For alpha = %1.2f, beta = %1.2f, gamma = %1.2f, RMSE is %3.4f MAPE is %3.2f" %(alpha_value, beta_value, gamma_value, rmse_opt, mape_opt))
For alpha = 0.88, beta = 0.00, gamma = 0.00, RMSE is 187.1334 MAPE is 32.89
tempResultsDf = pd.DataFrame({'Method': 'Holt_Winter', 'rmse': [rmse_opt], 'mape' : [mape_opt]})
resultsDf = pd.concat([resultsDf, tempResultsDf])
print(resultsDf)
Method rmse mape 0 RegressionOnTime 164.563000 27.94 0 RegressionOnTimeSeasonal 185.864000 32.59 0 Naive model 186.469000 31.72 0 Simple Average 279.196000 52.35 0 moving_avg_forecast_4 75.514000 13.05 0 moving_avg_forecast_6 88.576000 14.62 0 moving_avg_forecast_8 105.158000 17.03 0 moving_avg_forecast_12 134.071000 22.79 0 SES 203.706554 35.18 0 SES 203.706554 35.18 0 Holt_linear 164.871190 28.00 0 Holt_Winter 187.133365 32.89
y_hat_avg = test.copy()
model_fit = ExponentialSmoothing(np.asarray(train['Sales']) ,seasonal_periods = 12 ,trend='add', seasonal='Multiplicative').fit()
y_hat_avg['Holt_Winter_M'] = model_fit.forecast(len(test))
print('')
print('== Holt-Winters Additive ETS(A,A,M) Parameters ==')
print('')
alpha_value = np.round(pred.params['smoothing_level'], 4)
beta_value = np.round(model_fit.params['smoothing_trend'], 4)
gamma_value = np.round(model_fit.params['smoothing_seasonal'], 4)
print('Smoothing Level: ', alpha_value)
print('Smoothing Trend: ', beta_value)
print('Smoothing Seasonal: ', gamma_value)
print('Initial Level: ', np.round(model_fit.params['initial_level'], 4))
print('Initial Trend: ', np.round(model_fit.params['initial_trend'], 4))
print('Initial Seasons: ', np.round(model_fit.params['initial_seasons'], 4))
print('')
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\holtwinters\model.py:427: FutureWarning: After 0.13 initialization must be handled at model creation warnings.warn(
== Holt-Winters Additive ETS(A,A,M) Parameters == Smoothing Level: 0.8806 Smoothing Trend: 0.0 Smoothing Seasonal: 0.0 Initial Level: 180.789 Initial Trend: 6.2855 Initial Seasons: [0.9884 0.6238 0.8168 0.8347 0.7303 0.9035 0.8528 0.9835 0.885 1.0204 1.0109 0.9175]
plt.figure(figsize=(16,8))
plt.plot( train['Sales'], label='Train')
plt.plot(test['Sales'], label='Test')
plt.plot(y_hat_avg['Holt_Winter_M'], label='Holt_Winter_M')
plt.title('Holt-Winters Multiplicative Parameters:\n alpha = ' +
str(alpha_value) + ' Beta:' +
str(beta_value) +
' Gamma: ' + str(gamma_value))
plt.legend(loc='best')
<matplotlib.legend.Legend at 0x2126bba5fa0>
rmse_opt= np.sqrt(mean_squared_error(test['Sales'], y_hat_avg['Holt_Winter_M']))
mape_opt= MAPE(test['Sales'], y_hat_avg['Holt_Winter_M'])
print("For alpha = %1.2f, beta = %1.2f, gamma = %1.2f, RMSE is %3.4f MAPE is %3.2f" %(alpha_value, beta_value, gamma_value, rmse_opt, mape_opt))
For alpha = 0.88, beta = 0.00, gamma = 0.00, RMSE is 188.9831 MAPE is 33.03
tempResultsDf = pd.DataFrame({'Method': 'Holt_Winter M', 'rmse': [rmse_opt], 'mape' : [mape_opt]})
resultsDf = pd.concat([resultsDf, tempResultsDf])
print(resultsDf)
Method rmse mape 0 RegressionOnTime 164.563000 27.94 0 RegressionOnTimeSeasonal 185.864000 32.59 0 Naive model 186.469000 31.72 0 Simple Average 279.196000 52.35 0 moving_avg_forecast_4 75.514000 13.05 0 moving_avg_forecast_6 88.576000 14.62 0 moving_avg_forecast_8 105.158000 17.03 0 moving_avg_forecast_12 134.071000 22.79 0 SES 203.706554 35.18 0 SES 203.706554 35.18 0 Holt_linear 164.871190 28.00 0 Holt_Winter 187.133365 32.89 0 Holt_Winter M 188.983144 33.03
Window 4 크기의 이동 평균이 가장 성능이 좋았다...
자동 회귀( Auto-regressive, AR ), 이동 평균( Moving Average, MA ) 모델들은 예측에 자주 사용되는 모델들이다. AR, MA 모델을 합치면 ARMA라는 모델이 나온다. 초기의 ARMA와 ARIMA 모델은 1970년, Box and Jenkins가 개발했다.
ARMA 모델은 기본적으로 회귀 모델이다; 자동 회귀란 것은 여러 구간에서 변수들이 자체적으로 회귀를 돌려서 값이 나온다는 의미이다. AR 모델의 제 1가정은 시계열이 보합상태인 것이다. 안정화된 시계열 데이터라는 것은 mean, variance, autocorrelation 등의 지표가 모든 구간에서 동일하게 나온다는 의미이다.
시계열 데이터가 안정적이지 않으면, AR 모델 적용전에 데이터를 변형해야 한다.(???)
연속적인 관찰값과의 차이를 lag - 1 차이라고 불러보자. 계절성이 있는 시계열 데이터의 경우, Lag는 계절성 기간 동안에 나타날 것이다.
잔차의 White noise:
White noise는 잔차( $\epsilon_t$ )들의 결과이다. 서로 상관되어 있지 않고, 평균이 0이고, 표준편차가 상수인 정규 분포를 따른다. AR 모델의 제 1가정은 Error가 White noise를 따른다는 것이다.
Auto-Regression은 여러 시점에서 측정된 변수를 회귀하는 방법이다. Lag 1을 가지는 Auto-Regressive model은 AR(1)로 표현되고, 아래와 같은 성격을 가진다.
Augmented Dickey Fuller Test ( ADF )은 보합인지 확인하는 단위 Root 테스트이다.
서로 다른 Lag에 대한 자기 상관계수를 Plot 한 것을 ACF라 한다.
Plot은 관찰값들과 Lag 값들간의 상관관계를 함축해서 보여준다. X 축은 Lag, Y 축은 상관계수를 나타낸다. ( Plot 이 없는데...?, 작성하다가 누락한 듯.. )
서로 다른 Lag에 대한 부분 자기 상관 계수를 Plot 한것을 PACF라 한다.
관찰값과 Lag 값의 상관관계를 보여준다. Lag 값들은 이전의 Lag 값들에 사용되지 않는다.
두 그림( ??? 역시 그림 없음. ) 두 개 Plot은 95%, 99% 신뢰 구간에 따른 Bar chart 이다. The number of lags is p when:
* The model is AR if the ACF trails off after a lag and has a hard cut-off in the PACF after a lag. This lag is taken as the value for p.
* The model is MA if the PACF trails off after a lag and has a hard cut-off in the ACF after the lag. This lag value is taken as the value for q.
* The model is a mix of AR and MA if both the ACF and PACF trail off.
from statsmodels.tsa.stattools import adfuller
data= pd.read_csv(
os.path.join( DATA_ROOT, 'time-series-data', 'TractorSales.csv' ),
header=0, parse_dates=[0], squeeze=True
)
dates= pd.date_range(start='2003-01-01', freq='MS', periods=len(data))
data['Month']= dates.month
data['Month']= data['Month'].apply(lambda x: calendar.month_abbr[x])
data['Year']= dates.year
data.drop(['Month-Year'], axis=1, inplace=True)
data.rename(columns={'Number of Tractor Sold':'Tractor-Sales'}, inplace=True)
data= data[['Month', 'Year', 'Tractor-Sales']]
data.set_index(dates, inplace=True)
sales_ts = data['Tractor-Sales']
result = adfuller(sales_ts)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
ADF Statistic: 1.108825 p-value: 0.995291
ADF로 보아하니, 본 데이터는 Stationary 하지 않다. Differencing - Lag 1을 시도해 본다.
sales_ts_diff = sales_ts - sales_ts.shift(periods=1) # 1 뒤로 Lag
sales_ts_diff.dropna(inplace=True)
result = adfuller(sales_ts_diff)
pval = result[1]
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
if pval < 0.05:
print('Data is stationary')
else:
print('Data after differencing is not stationary; so try log diff')
sales_ts_log = np.log10(sales_ts)
sales_ts_log.dropna(inplace=True)
sales_ts_log_diff = sales_ts_log.diff(periods=1)
sales_ts_log_diff.dropna(inplace=True)
result = adfuller(sales_ts_log_diff)
pval = result[1]
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
if pval < 0.05:
print('Data after log differencing is stationary')
else:
print('Data after log differencing is not stationary; try second order differencing')
sales_ts_log_diff2 = sales_ts_log.diff(periods = 2)
sales_ts_log_diff2.dropna(inplace=True)
result = adfuller(sales_ts_log_diff2)
pval = result[1]
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
if pval < 0.05:
print('Data after log differencing 2nd order is stationary')
else:
print('Data after log differencing 2nd order is not stationary')
ADF Statistic: -2.543481 p-value: 0.105250 Data after differencing is not stationary; so try log diff ADF Statistic: -2.680467 p-value: 0.077480 Data after log differencing is not stationary; try second order differencing ADF Statistic: -3.200722 p-value: 0.019943 Data after log differencing 2nd order is stationary
2차례 Differencing 하고 나니 Stationary 해졌다....
#ACF and PACF plots:
from statsmodels.tsa.stattools import acf, pacf
import matplotlib.pyplot as plt
lag_acf = acf(sales_ts_log_diff2, nlags=20)
lag_pacf = pacf(sales_ts_log_diff2, nlags=20, method='ols')
#Plot ACF:
plt.figure(figsize = (15,5))
plt.subplot(121)
plt.stem(lag_acf)
plt.axhline(y = 0, linestyle='--',color='black')
plt.axhline(y = -1.96/np.sqrt(len(sales_ts_log_diff2)),linestyle='--',color='gray')
plt.axhline(y = 1.96/np.sqrt(len(sales_ts_log_diff2)),linestyle='--',color='gray')
plt.xticks(range(0,22,1))
plt.xlabel('Lag')
plt.ylabel('ACF')
plt.title('Autocorrelation Function')
#Plot PACF:
plt.subplot(122)
plt.stem(lag_pacf)
plt.axhline(y = 0, linestyle = '--', color = 'black')
plt.axhline(y =-1.96/np.sqrt(len(sales_ts_log_diff2)), linestyle = '--', color = 'gray')
plt.axhline(y = 1.96/np.sqrt(len(sales_ts_log_diff2)),linestyle = '--', color = 'gray')
plt.xlabel('Lag')
plt.xticks(range(0,22,1))
plt.ylabel('PACF')
plt.title('Partial Autocorrelation Function')
plt.tight_layout()
plt.show()
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\stattools.py:667: FutureWarning: fft=True will become the default after the release of the 0.12 release of statsmodels. To suppress this warning, explicitly set fft=False. warnings.warn(
Lag differeced 에 따른 ACF, PACF 결과를 살펴보자. ( Cut off의 의미를 알아야 하는데 잘 모르겠음..., 왜 Lag 5, 2가 선택되었는지 잘 모르겠음... )
d = 2로 설정한다. 2차례의 Differencing을 걸쳐 데이터가 Stationary 해졌기 때문이다.
random walk는 수학적 기법이다. Stochastic, Random process라고 알려져 있는데, 특정 Integer 구간에서 연속적인 Random steps로 구성된 것을 의미한다.
An elementary example of a random walk is the random walk on the integer number line, ${\displaystyle \mathbb {Z} } $, 0 에서 각 Step을 +1 이나 -1로 동일한 확률로 이동하는 것을 의미한다.
Random walk process is defined as: y(t) = $\beta_0 + \beta_1 X_{t-1}$ + $\epsilon_t$
where
이번 예제에서는, Random walk를 하기 위해서 동일한 Seed를 이용해 Random number를 만들어내는 함수를 이용할 것이다.
from random import seed, random
seed(1) # Seed 1 고정
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)
for i in range(1, 1000):
movement = -1 if random() < 0.5 else 1
value = random_walk[i-1] + movement
random_walk.append(value)
pd.plotting.autocorrelation_plot(random_walk)
plt.show()
### Check stationary property
result = adfuller(random_walk)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
print('\t%s: %.3f' % (key, value))
ADF Statistic: 0.341605 p-value: 0.979175 Critical Values: 1%: -3.437 5%: -2.864 10%: -2.568
# Random walk 예측을 위한 데이터셋 준비
training_size = int(len(random_walk) * 0.70)
training, test = random_walk[0 : training_size], random_walk[training_size:]
predictions = list()
hist = training[-1]
for i in range(len(test)):
yhat = hist
predictions.append(yhat)
hist = test[i]
rmse = np.sqrt(mean_squared_error(test, predictions))
print('\n\nPredicting a Random Walk \n RMSE: %.3f' % rmse)
Predicting a Random Walk RMSE: 1.000
주식 가격은 아래와 같이 시뮬레이션 한다.
Random walk with drift 시뮬레이션 결과를 그려보자.
# 생성, Seed는 1234 지정
seed(1234)
rw_steps = np.random.normal(loc = 0.001, scale = 0.01, size = 1000) + 1 # Mean, STDEV, Size
### 초기 값을 1로 설정
rw_steps[0] = 1
### 주식 가격 시뮬레이션
Price = rw_steps * np.cumprod(rw_steps)
Price = Price * 100
### 그리기
plt.plot(rw_steps)
plt.title("Simulated random walk with drift")
plt.show()
### Stationary 확인
result = adfuller(Price)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
print('\t%s: %.3f' % (key, value))
ADF Statistic: 2.422457 p-value: 0.999020 Critical Values: 1%: -3.437 5%: -2.864 10%: -2.568
### 예측
training_size = int(len(Price) * 0.70)
training, test = random_walk[0 : training_size], random_walk[training_size:]
predictions = list()
hist = training[-1]
for i in range(len(test)):
yhat = hist
predictions.append(yhat)
hist = test[i]
rmse = np.sqrt(mean_squared_error(test, predictions))
print('\n\nPredicting a Random Walk with drift \nRMSE: %.3f' % rmse)
Predicting a Random Walk with drift RMSE: 1.000
대표적으로 시계열 예측에 사용되는 방법 중 하나는 ARIMA 모델이다. Auto regressive + Moving Average 이다. AR, I, MA 로 분해해서 볼 수 있는데:
AR(p): Autoregressive model, p는 얼마나 Lagged된 시리즈를 예측할 것인지를 결정하는 파라미터
I(d): Differencing part, d는 Series를 Stationary하게 만들기 위해서 얼마나 많이 Differencing 해야 하는지에 대한 파라미터이다.
MA(q): Moving average model, q는 예측 수식에서 Lagged 예측 Error term의 개수이다.
SARIMA는 계절성이 들어간 ARIMA이고, 계절성이 있는 경우에 사용된다.
Let's see the example of data from quandl의 데이터를 살펴보자. 영란은행의 EURO - US Dollar Spot 환율 데이터이다.
from pandas import DataFrame
from io import StringIO
import time, json
from datetime import date
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.seasonal import seasonal_decompose
from sklearn.metrics import mean_squared_error
%matplotlib inline
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6
bank_df = pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'BOE-XUDLERD.csv' ) )
bank_df
| Date | Value | |
|---|---|---|
| 0 | 2017-11-09 | 0.8603 |
| 1 | 2017-11-08 | 0.8631 |
| 2 | 2017-11-07 | 0.8639 |
| 3 | 2017-11-06 | 0.8631 |
| 4 | 2017-11-03 | 0.8608 |
| ... | ... | ... |
| 10832 | 1975-01-08 | 0.7554 |
| 10833 | 1975-01-07 | 0.7510 |
| 10834 | 1975-01-06 | 0.7524 |
| 10835 | 1975-01-03 | 0.7585 |
| 10836 | 1975-01-02 | 0.7633 |
10837 rows × 2 columns
# 시간을 축으로 만듦
bank_df['Date'] = pd.to_datetime(bank_df['Date'])
indexed_df = bank_df.set_index('Date')
ts = indexed_df['Value']
ts.head(5)
Date 2017-11-09 0.8603 2017-11-08 0.8631 2017-11-07 0.8639 2017-11-06 0.8631 2017-11-03 0.8608 Name: Value, dtype: float64
#Visualize the raw data
plt.plot(ts)
[<matplotlib.lines.Line2D at 0x2126b17b490>]
# Resample the data as it contains too much variations
ts_week = ts.resample('W').mean() # 주 단위로 Resampling
plt.plot(ts_week)
[<matplotlib.lines.Line2D at 0x2126b20d130>]
def test_stationarity(timeseries):
#Determing rolling statistics
rolmean = timeseries.rolling(window=52,center=False).mean()
rolstd = timeseries.rolling(window=52,center=False).std()
#Plot rolling statistics:
orig = plt.plot(timeseries, color='blue',label='Original')
mean = plt.plot(rolmean, color='red', label='Rolling Mean')
std = plt.plot(rolstd, color='black', label = 'Rolling Std')
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
plt.show(block=False)
#Perform Dickey-Fuller test:
print('Results of Dickey-Fuller Test:')
dftest = adfuller(timeseries, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
print(dfoutput)
test_stationarity(ts_week)
Results of Dickey-Fuller Test: Test Statistic -2.076341 p-value 0.254134 #Lags Used 2.000000 Number of Observations Used 2234.000000 Critical Value (1%) -3.433281 Critical Value (5%) -2.862835 Critical Value (10%) -2.567459 dtype: float64
검정 통계량이 5% Critical value보다 크고, p-value는 0.05보다 크다. Moving average는 상수가 아이며, H0 가설은 기각할 수가 없다. 결국, 주간 데이터 역시 Stationary 하지 않다.
어떤 방법을 써서라도, Stationary로 만들어야 겠다...
Stationary 한 데이터로 만들기 위해서 Differencing을 수행한다. 수행 전에 Log 데이터로 변환하는 것이 낫겠다.
ts_week_log = np.log(ts_week)
ts_week_log_diff = ts_week_log - ts_week_log.shift()
plt.plot(ts_week_log_diff)
[<matplotlib.lines.Line2D at 0x2126bf58a60>]
# Again confirming with the dickey-fuller test
ts_week_log_diff.dropna(inplace=True)
test_stationarity(ts_week_log_diff)
Results of Dickey-Fuller Test: Test Statistic -36.590004 p-value 0.000000 #Lags Used 0.000000 Number of Observations Used 2235.000000 Critical Value (1%) -3.433279 Critical Value (5%) -2.862834 Critical Value (10%) -2.567459 dtype: float64
추론
로그로 만드니까, Critical value는 1%보다 작고, 99% 신뢰 구간에 들어온다. Stationary 한 모델에 쓸 수 있는 ARIMA를 사용할 수 있게 되었다.
#ACF and PACF
lag_acf = acf(ts_week_log_diff, nlags=10)
lag_pacf = pacf(ts_week_log_diff, nlags=10, method='ols')
#Plot ACF:
plt.subplot(121)
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-7.96/np.sqrt(len(ts_week_log_diff)),linestyle='--',color='gray')
plt.axhline(y=7.96/np.sqrt(len(ts_week_log_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function')
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\stattools.py:667: FutureWarning: fft=True will become the default after the release of the 0.12 release of statsmodels. To suppress this warning, explicitly set fft=False. warnings.warn(
Text(0.5, 1.0, 'Autocorrelation Function')
#Plot PACF:
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-7.96/np.sqrt(len(ts_week_log_diff)),linestyle='--',color='gray')
plt.axhline(y=7.96/np.sqrt(len(ts_week_log_diff)),linestyle='--',color='gray')
plt.title('Partial Autocorrelation Function')
plt.tight_layout()
ACF, PACF를 기반으로 추론.
p, q 값을 정해보자.
# Optimal values fot ARIMA(p,d,q) model are (2,1,1). Hence plot the ARIMA model using the value (2,1,1)
model = ARIMA(ts_week_log, order=(2, 1, 1))
results_ARIMA = model.fit(disp=-1)
plt.plot(ts_week_log_diff)
plt.plot(results_ARIMA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_ARIMA.fittedvalues-ts_week_log_diff)**2))
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
Text(0.5, 1.0, 'RSS: 0.2806')
# print the results of the ARIMA model and plot the residuals
print(results_ARIMA.summary())
# plot residual errors
residuals = DataFrame(results_ARIMA.resid)
residuals.plot(kind='kde')
print(residuals.describe())
ARIMA Model Results
==============================================================================
Dep. Variable: D.Value No. Observations: 2236
Model: ARIMA(2, 1, 1) Log Likelihood 6870.601
Method: css-mle S.D. of innovations 0.011
Date: Mon, 29 Nov 2021 AIC -13731.202
Time: 00:01:53 BIC -13702.640
Sample: 01-12-1975 HQIC -13720.773
- 11-12-2017
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
const 5.509e-05 0.000 0.178 0.859 -0.001 0.001
ar.L1.D.Value -0.0901 0.487 -0.185 0.853 -1.044 0.864
ar.L2.D.Value 0.0602 0.128 0.469 0.639 -0.191 0.312
ma.L1.D.Value 0.3476 0.485 0.716 0.474 -0.604 1.299
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 -3.3955 +0.0000j 3.3955 0.5000
AR.2 4.8928 +0.0000j 4.8928 0.0000
MA.1 -2.8770 +0.0000j 2.8770 0.5000
-----------------------------------------------------------------------------
0
count 2236.000000
mean 0.000001
std 0.011205
min -0.061306
25% -0.006725
50% -0.000228
75% 0.006869
max 0.064140
#Predictions
predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True)
print (predictions_ARIMA_diff.head())
Date 1975-01-12 0.000055 1975-01-19 -0.002420 1975-01-26 0.000987 1975-02-02 -0.004103 1975-02-09 -0.001134 Freq: W-SUN, dtype: float64
모델 우리가 원하는 결과를 보여주었다. Model의 예측 값을 원래 Scale로 복원해서 보자. 단, First order differencing을 제거하고, prediction을 원래 Scale로 만들기 위해서 Exponent 해야 한다.
predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum()
predictions_ARIMA_log = pd.Series(ts_week_log.iloc[0], index=ts_week_log.index)
predictions_ARIMA_log = predictions_ARIMA_log.add(predictions_ARIMA_diff_cumsum,fill_value=0)
predictions_ARIMA = np.exp(predictions_ARIMA_log)
plt.plot(ts_week)
plt.plot(predictions_ARIMA)
plt.title('RMSE: %.4f'% np.sqrt(sum((predictions_ARIMA-ts_week)**2)/len(ts_week)))
Text(0.5, 1.0, 'RMSE: 0.1353')
#Training and testing datsets
size = int(len(ts_week_log) - 15)
train, test = ts_week_log[0:size], ts_week_log[size:len(ts_week_log)]
history = [x for x in train]
predictions = list()
#Training the model and forecasting
size = int(len(ts_week_log) - 15)
train, test = ts_week_log[0:size], ts_week_log[size:len(ts_week_log)]
history = [x for x in train]
predictions = list()
print('Printing Predicted vs Expected Values...')
print('\n')
for t in range(len(test)):
model = ARIMA(history, order=(2,1,1))
model_fit = model.fit(disp=0)
output = model_fit.forecast()
yhat = output[0]
predictions.append(float(yhat))
obs = test[t]
history.append(obs)
print('predicted=%f, expected=%f' % (np.exp(yhat), np.exp(obs)))
Printing Predicted vs Expected Values...
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
predicted=0.860853, expected=0.862600
C:\Users\user\anaconda3\lib\site-packages\statsmodels\tsa\arima_model.py:472: FutureWarning:
statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have
been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .
between arima and model) and
statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.
statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and
is both well tested and maintained.
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
#Validating the model
error = mean_squared_error(test, predictions)
print('\n')
print('Printing Mean Squared Error of Predictions...')
print('Test MSE: %.6f' % error)
predictions_series = pd.Series(predictions, index = test.index)
Printing Mean Squared Error of Predictions... Test MSE: 0.000043
#Plotting forecasted vs Observed values
fig, ax = plt.subplots()
ax.set(title='Spot Exchange Rate, Euro into USD', xlabel='Date', ylabel='Euro into USD')
ax.plot(ts_week[-60:], 'o', label='observed')
ax.plot(np.exp(predictions_series), 'g', label='rolling one-step out-of-sample forecast')
legend = ax.legend(loc='upper left')
legend.get_frame().set_facecolor('w')
# 오.... 꽤나 잘 맞음....
이것은 서로 다른 조합의 p, q, d 값들을 이용하고, 이의 최적값을 만들기 위해 AIC (Akaike Information Criterion)와 BIC(Bayesian Information Criterion) values 값을 가용한다.
ARIMA 모델에서 사용한 3개 파라미터 p, q, d는 모델 입장에서 Seasonaliry, Trend, Noise에 대한 중요한 Aspect가 된다.
Seasonal ARIMA 모델을 사용한다면, P, Q, D는 모델의 계절성 컴포넌트와 관계된다.
auto.arima 모듈은 단일 시계열 데이터에 대해서 AIC, AICc, BIC 값을 이용해서 최적읜 Model을 만들어 낸다. 이 함수는 주어진 제약사항 내에서 가능한 모델을 전부 검색하는 기능을 수행한다.
The Akaike information criterion (AIC)는 주어진 데이터에 대한 통계 모델의 상대적인 품질을 측정하는 방법이다. AIC를 통해서 모델을 선택할 수 있다.
Bayesian information criterion (BIC) 혹은 Schwarz information criterion (also SIC, SBC, SBIC)은 유한한 세트의 모델에서 모델을 선택하는 어떤 기준이다. BIC 값이 가장 작은 것을 선택하는 것이 좋다. Likelihood function이 사용되는데 AIC와 연관되어 있다.
트랙터 판매량 데이터를 auto.arima로 돌려보자.
import sys, warnings, itertools
warnings.filterwarnings("ignore")
import statsmodels.api as sm, \
statsmodels.tsa.api as sm, \
statsmodels.formula.api as smf
tractor_sales_Series = pd.read_csv( os.path.join( DATA_ROOT, 'time-series-data', 'TractorSales.csv' ) )
tractor_sales_Series.head(5)
| Month-Year | Number of Tractor Sold | |
|---|---|---|
| 0 | 3-Jan | 141 |
| 1 | 3-Feb | 157 |
| 2 | 3-Mar | 185 |
| 3 | 3-Apr | 199 |
| 4 | 3-May | 203 |
dates = pd.date_range(start='2003-01-01', freq='MS', periods=len(tractor_sales_Series))
import calendar
data['Month'] = dates.month
data['Month'] = data['Month'].apply(lambda x: calendar.month_abbr[x])
data['Year'] = dates.year
#data.drop(['Month-Year'], axis=1, inplace=True)
data.rename(columns={'Number of Tractor Sold':'Tractor-Sales'}, inplace=True)
data = data[['Month', 'Year', 'Tractor-Sales']]
data.set_index(dates, inplace=True)
data.head(5)
| Month | Year | Tractor-Sales | |
|---|---|---|---|
| 2003-01-01 | Jan | 2003 | 141 |
| 2003-02-01 | Feb | 2003 | 157 |
| 2003-03-01 | Mar | 2003 | 185 |
| 2003-04-01 | Apr | 2003 | 199 |
| 2003-05-01 | May | 2003 | 203 |
# extract out the time-series
sales_ts = data['Tractor-Sales']
plt.figure(figsize=(8, 4))
plt.plot(sales_ts)
plt.xlabel('Years')
plt.ylabel('Tractor Sales')
Text(0, 0.5, 'Tractor Sales')
추세와 Multiplicative 계절성이 있는 것이 보여진다. Window size 4, 6, 8, 12로 MA를 돌려보자.
fig, axes = plt.subplots(2, 2, sharey=False, sharex=False)
fig.set_figwidth(14)
fig.set_figheight(8)
axes[0][0].plot(sales_ts.index, sales_ts, label='Original')
axes[0][0].plot(sales_ts.index, sales_ts.rolling(window=4).mean(), label='4-Months Rolling Mean')
axes[0][0].set_xlabel("Years")
axes[0][0].set_ylabel("Number of Tractor's Sold")
axes[0][0].set_title("4-Months Moving Average")
axes[0][0].legend(loc='best')
axes[0][1].plot(sales_ts.index, sales_ts, label='Original')
axes[0][1].plot(sales_ts.index, sales_ts.rolling(window=6).mean(), label='6-Months Rolling Mean')
axes[0][1].set_xlabel("Years")
axes[0][1].set_ylabel("Number of Tractor's Sold")
axes[0][1].set_title("6-Months Moving Average")
axes[0][1].legend(loc='best')
axes[1][0].plot(sales_ts.index, sales_ts, label='Original')
axes[1][0].plot(sales_ts.index, sales_ts.rolling(window=8).mean(), label='8-Months Rolling Mean')
axes[1][0].set_xlabel("Years")
axes[1][0].set_ylabel("Number of Tractor's Sold")
axes[1][0].set_title("8-Months Moving Average")
axes[1][0].legend(loc='best')
axes[1][1].plot(sales_ts.index, sales_ts, label='Original')
axes[1][1].plot(sales_ts.index, sales_ts.rolling(window=12).mean(), label='12-Months Rolling Mean')
axes[1][1].set_xlabel("Years")
axes[1][1].set_ylabel("Number of Tractor's Sold")
axes[1][1].set_title("12-Months Moving Average")
axes[1][1].legend(loc='best')
plt.tight_layout()
plt.show()
#Determing rolling statistics
rolmean = sales_ts.rolling(window = 4).mean()
rolstd = sales_ts.rolling(window = 4).std()
#Plot rolling statistics:
orig = plt.plot(sales_ts, label='Original')
mean = plt.plot(rolmean, label='Rolling Mean')
std = plt.plot(rolstd, label = 'Rolling Std')
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
Text(0.5, 1.0, 'Rolling Mean & Standard Deviation')
Dickey-Fuller Test - Let's run the Dicky Fuller Test on the timeseries and verify the null hypothesis that the TS is non-stationary.
from statsmodels.tsa.stattools import adfuller
dftest = adfuller(sales_ts)
dftest
print('DF test statistic is %3.3f' %dftest[0])
print('DF test p-value is %1.4f' %dftest[1])
DF test statistic is 1.109 DF test p-value is 0.9953
표준편차의 분산이 작지만, Rolling mean이 시간에 따라 증가하므로, Stationary 하지 않다. 검정 통계량 역시 Critical value 보다 높다. 위 그래프로 보아 월별 패턴이 있는 것으로 보여진다. Seasonal component를 분해해보자.
월별로 트랙터 판매량이 얼마나 변화했는지를 보자. 연도별 누적 Plot으로 확인해 볼 수 있다.
monthly_sales_data = pd.pivot_table(data, values = "Tractor-Sales", columns = "Year", index = "Month")
monthly_sales_data
| Year | 2003 | 2004 | 2005 | 2006 | 2007 | 2008 | 2009 | 2010 | 2011 | 2012 | 2013 | 2014 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Month | ||||||||||||
| Apr | 199 | 208 | 251 | 279 | 362 | 350 | 414 | 482 | 536 | 536 | 610 | 710 |
| Aug | 207 | 238 | 279 | 339 | 381 | 410 | 486 | 567 | 654 | 707 | 783 | 848 |
| Dec | 165 | 196 | 232 | 272 | 281 | 321 | 389 | 428 | 470 | 472 | 567 | 605 |
| Feb | 157 | 168 | 200 | 239 | 261 | 250 | 310 | 368 | 400 | 423 | 455 | 520 |
| Jan | 141 | 145 | 183 | 215 | 247 | 257 | 305 | 358 | 397 | 428 | 454 | 525 |
| Jul | 207 | 238 | 279 | 322 | 370 | 423 | 510 | 578 | 651 | 687 | 767 | 871 |
| Jun | 189 | 209 | 249 | 305 | 340 | 370 | 441 | 524 | 591 | 609 | 661 | 749 |
| Mar | 185 | 197 | 249 | 270 | 330 | 329 | 374 | 444 | 498 | 507 | 568 | 587 |
| May | 203 | 210 | 289 | 307 | 385 | 393 | 454 | 534 | 596 | 610 | 706 | 793 |
| Nov | 138 | 152 | 194 | 229 | 239 | 270 | 315 | 360 | 406 | 412 | 481 | 519 |
| Oct | 150 | 168 | 204 | 241 | 266 | 289 | 345 | 386 | 437 | 452 | 513 | 581 |
| Sep | 171 | 199 | 232 | 263 | 299 | 326 | 393 | 447 | 509 | 509 | 583 | 640 |
monthly_sales_data = monthly_sales_data.reindex(index = ['Jan','Feb','Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])
monthly_sales_data
| Year | 2003 | 2004 | 2005 | 2006 | 2007 | 2008 | 2009 | 2010 | 2011 | 2012 | 2013 | 2014 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Month | ||||||||||||
| Jan | 141 | 145 | 183 | 215 | 247 | 257 | 305 | 358 | 397 | 428 | 454 | 525 |
| Feb | 157 | 168 | 200 | 239 | 261 | 250 | 310 | 368 | 400 | 423 | 455 | 520 |
| Mar | 185 | 197 | 249 | 270 | 330 | 329 | 374 | 444 | 498 | 507 | 568 | 587 |
| Apr | 199 | 208 | 251 | 279 | 362 | 350 | 414 | 482 | 536 | 536 | 610 | 710 |
| May | 203 | 210 | 289 | 307 | 385 | 393 | 454 | 534 | 596 | 610 | 706 | 793 |
| Jun | 189 | 209 | 249 | 305 | 340 | 370 | 441 | 524 | 591 | 609 | 661 | 749 |
| Jul | 207 | 238 | 279 | 322 | 370 | 423 | 510 | 578 | 651 | 687 | 767 | 871 |
| Aug | 207 | 238 | 279 | 339 | 381 | 410 | 486 | 567 | 654 | 707 | 783 | 848 |
| Sep | 171 | 199 | 232 | 263 | 299 | 326 | 393 | 447 | 509 | 509 | 583 | 640 |
| Oct | 150 | 168 | 204 | 241 | 266 | 289 | 345 | 386 | 437 | 452 | 513 | 581 |
| Nov | 138 | 152 | 194 | 229 | 239 | 270 | 315 | 360 | 406 | 412 | 481 | 519 |
| Dec | 165 | 196 | 232 | 272 | 281 | 321 | 389 | 428 | 470 | 472 | 567 | 605 |
monthly_sales_data.plot()
<AxesSubplot:xlabel='Month'>
yearly_sales_data = pd.pivot_table(data, values = "Tractor-Sales", columns = "Month", index = "Year")
yearly_sales_data = yearly_sales_data[['Jan','Feb','Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']]
yearly_sales_data
| Month | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Year | ||||||||||||
| 2003 | 141 | 157 | 185 | 199 | 203 | 189 | 207 | 207 | 171 | 150 | 138 | 165 |
| 2004 | 145 | 168 | 197 | 208 | 210 | 209 | 238 | 238 | 199 | 168 | 152 | 196 |
| 2005 | 183 | 200 | 249 | 251 | 289 | 249 | 279 | 279 | 232 | 204 | 194 | 232 |
| 2006 | 215 | 239 | 270 | 279 | 307 | 305 | 322 | 339 | 263 | 241 | 229 | 272 |
| 2007 | 247 | 261 | 330 | 362 | 385 | 340 | 370 | 381 | 299 | 266 | 239 | 281 |
| 2008 | 257 | 250 | 329 | 350 | 393 | 370 | 423 | 410 | 326 | 289 | 270 | 321 |
| 2009 | 305 | 310 | 374 | 414 | 454 | 441 | 510 | 486 | 393 | 345 | 315 | 389 |
| 2010 | 358 | 368 | 444 | 482 | 534 | 524 | 578 | 567 | 447 | 386 | 360 | 428 |
| 2011 | 397 | 400 | 498 | 536 | 596 | 591 | 651 | 654 | 509 | 437 | 406 | 470 |
| 2012 | 428 | 423 | 507 | 536 | 610 | 609 | 687 | 707 | 509 | 452 | 412 | 472 |
| 2013 | 454 | 455 | 568 | 610 | 706 | 661 | 767 | 783 | 583 | 513 | 481 | 567 |
| 2014 | 525 | 520 | 587 | 710 | 793 | 749 | 871 | 848 | 640 | 581 | 519 | 605 |
yearly_sales_data.plot()
<AxesSubplot:xlabel='Year'>
yearly_sales_data.boxplot()
<AxesSubplot:>
트랙터 판매량은 매년 감소 없이 꾸준히 증가한다. 7, 8월 판매량이 Peak 이고, 해당 월의 변동폭도 다른 월보다 크다. 월별 Mean 값을 보면 1년 동안 판매량이 늘었다가 줄어드는 추세가 보여진다. Seasonal effect가 12개월에 걸친 Cycle이란 얘기이다.
decomposition = sm.tsa.seasonal.seasonal_decompose(sales_ts, model='multiplicative')
fig = decomposition.plot()
fig.set_figwidth(8)
fig.set_figheight(6)
fig.suptitle('Decomposition of multiplicative time series')
plt.show()
1) Trend( 추세 ): 12개월 MA가 직선과 유사해서, 선형 회귀 방법을 사용할 수 있다.
2) Seasonality( 계절성 ): 월별로 일정한 패턴을 보인다. 월별 Seasonal component는 월별 Trend 값을 뺀 값이다. 추세 값은 아래 수식으로 제거된다.
Seasonality_t × Remainder_t = Y_t/Trend_t
3) Irregular Remainder (random): 추세와 계절성을 제거한 후 잔차 값이다. 계산은 아래 식으로 이루어진다.
Remainder_t = Y_t / (Trend_t × Seasonality_t)
plt.figure(figsize=(8, 4))
plt.plot(sales_ts.diff(periods=1))
plt.xlabel('Years')
plt.ylabel('Tractor Sales')
Text(0, 0.5, 'Tractor Sales')
We observe seasonality even after differencing.
plt.figure(figsize=(8, 4))
plt.plot(np.log10(sales_ts))
plt.xlabel('Years')
plt.ylabel('Log (Tractor Sales)')
Text(0, 0.5, 'Log (Tractor Sales)')
We observe trend and seasonality even after taking log of the observations.
plt.figure(figsize=(10, 5))
plt.plot(np.log10(sales_ts).diff(periods=1))
plt.xlabel('Years')
plt.ylabel('Differenced Log (Tractor Sales)')
Text(0, 0.5, 'Differenced Log (Tractor Sales)')
sales_ts_log = np.log10(sales_ts)
sales_ts_log.dropna(inplace=True)
sales_ts_log_diff = sales_ts_log.diff(periods=1) # same as ts_log_diff = ts_log - ts_log.shift(periods=1)
sales_ts_log_diff.dropna(inplace=True)
fig, axes = plt.subplots(1, 2)
fig.set_figwidth(12)
fig.set_figheight(4)
sm.graphics.tsaplots.plot_acf(sales_ts_log, lags=30, ax=axes[0])
sm.graphics.tsaplots.plot_pacf(sales_ts_log, lags=30, ax=axes[1])
plt.tight_layout()
Nonstationary series have an ACF that remains significant for half a dozen or more lags, rather than quickly declining to zero. You must difference such a series until it is stationary before you can identify the process
The above ACF is “decaying”, or decreasing, very slowly, and remains well above the significance range (blue band) for at least a dozen lags. This is indicative of a non-stationary series.
fig, axes = plt.subplots(1, 2)
fig.set_figwidth(12)
fig.set_figheight(4)
plt.xticks(range(0,30,1), rotation = 90)
sm.graphics.tsaplots.plot_acf(sales_ts_log_diff, lags=30, ax=axes[0])
sm.graphics.tsaplots.plot_pacf(sales_ts_log_diff, lags=30, ax=axes[1])
plt.tight_layout()
The above ACF has “decayed” fast and remains within the significance range (blue band) except for a few (5) lags. This is indicative of a stationary series.
# Define the p, d and q parameters to take any value between 0 and 2
p = d = q = range(0, 2)
# Generate all different combinations of p, d and q triplets
pdq = list(itertools.product(p, d, q))
# Generate all different combinations of seasonal p, q and q triplets
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
pdq
[(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1)]
seasonal_pdq
[(0, 0, 0, 12), (0, 0, 1, 12), (0, 1, 0, 12), (0, 1, 1, 12), (1, 0, 0, 12), (1, 0, 1, 12), (1, 1, 0, 12), (1, 1, 1, 12)]
#Separate data into train and test
data['date'] = data.index
train = data[data.index < '2013-01-01']
test = data[data.index >= '2013-01-01']
train_sales_ts_log = np.log10(train['Tractor-Sales'])
best_aic = np.inf
best_pdq = None
best_seasonal_pdq = None
temp_model = None
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
temp_model = sm.tsa.statespace.SARIMAX(train_sales_ts_log,
order = param,
seasonal_order = param_seasonal,
enforce_stationarity=True)
results = temp_model.fit()
if results.aic < best_aic:
best_aic = results.aic
best_pdq = param
best_seasonal_pdq = param_seasonal
except:
#print("Unexpected error:", sys.exc_info()[0])
continue
print("Best SARIMAX{}x{}12 model - AIC:{}".format(best_pdq, best_seasonal_pdq, best_aic))
Best SARIMAXNonexNone12 model - AIC:inf
For ARIMA(p, d, q) × (P, D, Q)S, we got SARIMAX(0, 1, 1)x(1, 0, 1, 12)12 model with the least AIC:-600.0908420381976
Here,
best_model = sm.tsa.statespace.sarimax.SARIMAX(train_sales_ts_log,
order=(0, 1, 1),
seasonal_order=(1, 0, 1, 12),
enforce_stationarity=True)
best_results = best_model.fit()
print(best_results.summary().tables[0])
print(best_results.summary().tables[1])
SARIMAX Results
==========================================================================================
Dep. Variable: Tractor-Sales No. Observations: 120
Model: SARIMAX(0, 1, 1)x(1, 0, 1, 12) Log Likelihood 304.045
Date: Mon, 29 Nov 2021 AIC -600.091
Time: 00:39:13 BIC -588.974
Sample: 01-01-2003 HQIC -595.577
- 12-01-2012
Covariance Type: opg
==========================================================================================
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ma.L1 -0.2954 0.080 -3.715 0.000 -0.451 -0.140
ar.S.L12 0.9919 0.007 142.324 0.000 0.978 1.006
ma.S.L12 -0.5397 0.105 -5.120 0.000 -0.746 -0.333
sigma2 0.0003 3.29e-05 7.979 0.000 0.000 0.000
==============================================================================
pred_dynamic = best_results.get_prediction(start=pd.to_datetime('2012-01-01'), dynamic=True, full_results=True)
pred_dynamic_ci = pred_dynamic.conf_int()
pred99 = best_results.get_forecast(steps=24, alpha=0.1)
# Extract the predicted and true values of our time series
sales_ts_forecasted = pred_dynamic.predicted_mean
testCopy = test.copy()
testCopy['sales_ts_forecasted'] = np.power(10, pred99.predicted_mean)
testCopy
| Month | Year | Tractor-Sales | date | sales_ts_forecasted | |
|---|---|---|---|---|---|
| 2013-01-01 | Jan | 2013 | 454 | 2013-01-01 | 438.763155 |
| 2013-02-01 | Feb | 2013 | 455 | 2013-02-01 | 440.518398 |
| 2013-03-01 | Mar | 2013 | 568 | 2013-03-01 | 534.141272 |
| 2013-04-01 | Apr | 2013 | 610 | 2013-04-01 | 570.737746 |
| 2013-05-01 | May | 2013 | 706 | 2013-05-01 | 638.525991 |
| 2013-06-01 | Jun | 2013 | 661 | 2013-06-01 | 629.978456 |
| 2013-07-01 | Jul | 2013 | 767 | 2013-07-01 | 703.459574 |
| 2013-08-01 | Aug | 2013 | 783 | 2013-08-01 | 709.204309 |
| 2013-09-01 | Sep | 2013 | 583 | 2013-09-01 | 538.542543 |
| 2013-10-01 | Oct | 2013 | 513 | 2013-10-01 | 473.187733 |
| 2013-11-01 | Nov | 2013 | 481 | 2013-11-01 | 435.999552 |
| 2013-12-01 | Dec | 2013 | 567 | 2013-12-01 | 506.322176 |
| 2014-01-01 | Jan | 2014 | 525 | 2014-01-01 | 468.481833 |
| 2014-02-01 | Feb | 2014 | 520 | 2014-02-01 | 470.340826 |
| 2014-03-01 | Mar | 2014 | 587 | 2014-03-01 | 569.416570 |
| 2014-04-01 | Apr | 2014 | 710 | 2014-04-01 | 608.104971 |
| 2014-05-01 | May | 2014 | 793 | 2014-05-01 | 679.716194 |
| 2014-06-01 | Jun | 2014 | 749 | 2014-06-01 | 670.690132 |
| 2014-07-01 | Jul | 2014 | 871 | 2014-07-01 | 748.254129 |
| 2014-08-01 | Aug | 2014 | 848 | 2014-08-01 | 754.315216 |
| 2014-09-01 | Sep | 2014 | 640 | 2014-09-01 | 574.070529 |
| 2014-10-01 | Oct | 2014 | 581 | 2014-10-01 | 504.930548 |
| 2014-11-01 | Nov | 2014 | 519 | 2014-11-01 | 465.554756 |
| 2014-12-01 | Dec | 2014 | 605 | 2014-12-01 | 539.993050 |
# Compute the root mean square error
mse = ((testCopy['Tractor-Sales'] - testCopy['sales_ts_forecasted']) ** 2).mean()
rmse = np.sqrt(mse)
print('The Root Mean Squared Error of our forecasts is {}'.format(round(rmse, 3)))
The Root Mean Squared Error of our forecasts is 65.756
axis = train['Tractor-Sales'].plot(label='Train Sales', figsize=(10, 6))
testCopy['Tractor-Sales'].plot(ax=axis, label='Test Sales', alpha=0.7)
testCopy['sales_ts_forecasted'].plot(ax=axis, label='Forecasted Sales', alpha=0.7)
axis.set_xlabel('Years')
axis.set_ylabel('Tractor Sales')
plt.legend(loc='best')
plt.show()
plt.close()
# Get forecast 36 steps (3 years) ahead in future
n_steps = 36
pred_uc_99 = best_results.get_forecast(steps=36, alpha=0.01) # alpha=0.01 signifies 99% confidence interval
pred_uc_95 = best_results.get_forecast(steps=36, alpha=0.05) # alpha=0.05 95% CI
# Get confidence intervals 95% & 99% of the forecasts
pred_ci_99 = pred_uc_99.conf_int()
pred_ci_95 = pred_uc_95.conf_int()
n_steps = 36
idx = pd.date_range(data.index[-1], periods=n_steps, freq='MS')
fc_95 = pd.DataFrame(np.column_stack([np.power(10, pred_uc_95.predicted_mean), np.power(10, pred_ci_95)]),
index=idx, columns=['forecast', 'lower_ci_95', 'upper_ci_95'])
fc_99 = pd.DataFrame(np.column_stack([np.power(10, pred_ci_99)]),
index=idx, columns=['lower_ci_99', 'upper_ci_99'])
fc_all = fc_95.combine_first(fc_99)
fc_all = fc_all[['forecast', 'lower_ci_95', 'upper_ci_95', 'lower_ci_99', 'upper_ci_99']] # just reordering columns
fc_all.head()
| forecast | lower_ci_95 | upper_ci_95 | lower_ci_99 | upper_ci_99 | |
|---|---|---|---|---|---|
| 2014-12-01 | 438.763155 | 407.819738 | 472.054410 | 407.819738 | 472.054410 |
| 2015-01-01 | 440.518398 | 402.819122 | 481.745897 | 402.819122 | 481.745897 |
| 2015-02-01 | 534.141272 | 481.746054 | 592.235050 | 481.746054 | 592.235050 |
| 2015-03-01 | 570.737746 | 508.538881 | 640.544090 | 508.538881 | 640.544090 |
| 2015-04-01 | 638.525991 | 562.725074 | 724.537541 | 562.725074 | 724.537541 |
# plot the forecast along with the confidence band
axis = sales_ts.plot(label='Observed', figsize=(8, 4))
fc_all['forecast'].plot(ax=axis, label='Forecast', alpha=0.7)
axis.fill_between(fc_all.index, fc_all['lower_ci_95'], fc_all['upper_ci_95'], color='k', alpha=.15)
axis.set_xlabel('Years')
axis.set_ylabel('Tractor Sales')
plt.legend(loc='best')
plt.show()
best_results.plot_diagnostics(lags=30, figsize=(16,12))
plt.show()
We need to ensure that the residuals of our model are uncorrelated and normally distributed with zero-mean. If it is not that it signifies that the model can be further improved and we repeat the process with the residuals.
In this case, our model diagnostics suggests that the model residuals are normally distributed based on the following:
Those observations coupled with the fact that there are no spikes outside the insignificant zone for both ACF and PACF plots lead us to conclude that that residuals are random with no information or juice in them and our model produces a satisfactory fit that could help us understand our time series data and forecast future values. It sems that our ARIMA model is working fine.